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FORTRAN  SUBROUTINES  TO  SOLVE  THE  LINEAR  LEAST-SQUARES 
PROBLEM  AND  COMPUTE  THE  COMPLETE  ORTHOGONAL  FACTORIZATION 

by 

Margaret  H.  Wright  and  Steven  C.  Classman 

1 . Int  roductlon 
I . I . Overview 

The  topics  of  Interest  In  this  report  are: 

(1)  The  linear  least-squares  problem  — given  a real  m by  n 

matrix  A and  an  m-vector  b,  find  an  n-vcctor  x such  that 
2 

li  Ax  - b II 2 Is  a minimum;  this  problem  Includes  the  solution  of 
non-singular,  over-  and  under-determined  linear  systems; 

(11)  computation  of  the  complete  orthogonal  factorization  of  a real 
matrix  A of  rank  r — find  orthogonal  matrices  Q and  V, 
and  a non-singular  r by  r upper  triangular  matrix  R, 
such  that: 


Tlie  software  to  be  described  is  a modular  set  of  portable 


I 


Fortran  subroutines  designed  to  carry  out  the  computations  neces- 
sary to  solve  (i)  and  (il),  and  also  to  perform  several  related 
tasks.  These  routines  have  already  been  used  in  several  applica- 
tions — e.g.,  in  large-scale  linear  programming  (Perold  and  Dantzig, 
1978),  in  solving  ill-conditioned  linear  systems  arising  from 
stochastic  processes,  and  in  qtiadratic  programming. 

If  tlie  reader  is  interested  only  in  using  the  software,  the 
descriptions  and  documentation  are  contained  in  Section  6. 

1.2.  Householde r Transformations 

The  techniques  used  to  solve  (1)  and  (il)  are  based  on  the 
construction  of  a sequence  of  orthogonal  transformations  designed 
to  reduce  the  original  matrix  to  a special  form.  The  advantages  of 
this  approach  in  terms  of  elegance,  efficiency  and  numerical  stabil- 
ity are  well  known.  The  theory  and  procedures  were  popularized  pri- 
marily by  G.H.  Golub,  and  are  explained  in  detail  in  numerous  refer- 
ences (e.g.,  Lawson  and  Hanson,  1974;  Stewart,  1973);  only  a brief 
description  will  be  given  here. 

For  any  non-zero  vector  u,  we  define  the  corresponding 
Housetu5l^e_r  trans formation  (or  Hou£ehoH^r  m^tr^)  as  an  elementary 
matrix  of  t!ie  form: 


T T 

H(u)  F I - I - ^ 

llul|2 


2 


In  this  context,  the  vector  u Is  called  the  Householder  vector 


and  the  scalar  P is  termed  the  scaling  factor  for  the  transforma- 
tion . 

Householder  matrices  are  symmetric  and  orthogonal  (so  that 
Euclidean  length  is  preserved  by  their  application).  In  addition, 
tliey  have  the  following  useful  properties: 

(a)  for  any  two  distinct  vectors  a and  b of  equal  Euclidean 
length,  there  exists  a Householder  matrix  that  will  transform 
one  into  the  other.  In  order  to  construct  such  a transforma- 
tion H,  we  seek  a vector  u that  satisfies: 

T 

Ha  = (I  - ■'^)a  = b , 

so  that 

T 

(-  ^)u  = b - a . (1) 

From  (1),  u must  be  a vector  in  the  direction  (b  - a) , and 
is  non-zero  because  a and  b are  distinct. 

(b)  The  vector  that  results  from  applying  a Householder  transforma- 
tion is  of  a special  form,  since  for  any  vector  c: 

He  = (I  - ^)c  = c - u(^)  . (2) 
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Tile  transformed  vector  He  Is  thus  given  by  the  difference 

between  the  original  vector  and  a multiple  of  the  Householder 

vector.  Consequently,  the  transformed  vector  is  identical 

to  the  original  vector  in  all  components  where  the  Householder 

vector  is  zero;  furthermore,  the  vector  c is  not  altered  at 

all  by  the  Householder  transformation  if  it  is  orthogonal  to 

T 

the  Householder  vector,  i.e.,  the  inner  product  u c is  zero. 
Tliese  properties  of  Houseliolder  transformations  can  be  exploited 
to  reduce  a general  real  matrix  to  a form  that  permits  easy  solution 
of  problems  (1)  and  (li). 


U 


♦ 

2 . The  F ul 1-Rank  Linear  Least-Squares  Problem 

In  this  section,  it  will  be  assumed  that  rank (A;  = n,  so 
that  m n.  We  defer  until  Section  A the  complex  issue  of  deter- 
mining the  rank . 

2.1.  Reduction  to  Upper  Triangular  Form  by  Transformat ion  From 
the  Left 

Because  of  properties  (a)  and  (b),  it  is  possible  to  construct 
a sequence  of  Householder  matrices  — H^,  ...,  — to  be  applied 

to  the  original  matrix  on  the  left  to  reduce  it  to  upper  triangular 
fqrnn,  i .e. , 

“nVl  W'’* 

where  R is  an  n by  n upper  triangular  matrix,  and  Q is  an 
m by  m orthogonal  matrix  (the  product  of  the  matrices  •••  H^^). 

The  fir.st  step  of  this  reduction  is  to  construct  a Householder  | 

i 

transformation  to  annihilate  components  2 through  m of  the  ] 

first  column  of  A.  From  (1),  it  follows  that  u^^ , the  vector  cor- 
responding to  Hj^,  will  be  equal  to  the  first  column  of  A,  except 
for  the  first  component.  Because  Householder  transformations  preserve 

F.uclidean  length,  the  transformed  first  column  will  be  given  by 

T 22  21/2 

0 •••  0]  , where  + ...  + . Using  (1), 


the  vector  Uj^  is  given  by: 


r 


1 


n 


^11  “ *^11 
"*21 

’*1  = . 

/nd 

To  avoid  cancellation  error  in  computing  u^,  the  sign  of  r^^^^  is 
chosen  to  be  opposite  that  of  the  original  thus,  the  first 

component  of  u^  is  "*■  sign (a^^)  I r^j^  1 ) . 

After  application  of  , the  first  column  of  the  transformed 
matrix  has  the  form  of  the  first  column  of  an  upper  triangular  matrix. 
For  example,  if  m = 4,  n = 3,  the  transformed  matrix  is  given  by: 


’^ll 

^12 

"*13 

0 

^22 

^23 

0 

^32 

^33 

0 

^42 

^43 

and,  in  general,  all  elements  of  A are  altered. 

In  constructing  the  second  Householder  transformation  H2, 
the  aim  is  to  annihilate  components  3 through  m of  the  second 
column,  while  leaving  the  first  column  unaltered.  From  (2),  this  can 
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I 

r 


be  accomplished  by  setting  the  first  component  of  the  Householder 
vector  to  zero,  for  ttien  is  ortliogonal  to  the  first  column 

of  H|A.  Furthermore,  because  the  first  component  of  u^  is  zero, 
the  first  row  of  H^A  will  not  be  altered  by  H2.  Since  the  first 
row  and  column  remain  unchanged,  they  can  be  ignored;  in  essence, 
the  second  transformation  is  applied  to  the  "remaining  matrix"  of 
(m  - 1)  rows  and  (n  - 1)  columns  that  results  from  omitting  the  first 
row  and  column  of  Hj^A. 

After  H,  is  applied  to  H^^A,  the  second  column  is  in  the 
desired  form.  With  the  example  above, 


’'11 

^2 

^13 

0 

r22 

^23 

0 

0 

"^33 

0 

0 

"^43 

where  r22  tind  the  doubly  barred  elements  have  been  altered  by  H2. 

Tills  process  is  continued  until  the  matrix  lias  been  reduced  to 
upper  triangular  form.  Each  step  may  be  viewed  as  transforming  the 
"first"  colum.n  of  a successively  smaller  remaining  matrix,  since  by 
construction  the  (k  + l)-st  transformation  does  not  alter  rows  or 
columns  1 through  k.  It  should  be  noted  that  the  assumption  of 
full  rank  is  essential  to  ensure  that  the  first  column  of  the  remain- 
ing matrix  Is  never  identically  zero  at  any  stage  of  the  reduction. 


r 
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Details  of  the  actual  computation,  operation  counts,  and 
data  structures  used  in  cc.  structing,  applying  and  storing  these 
transformations  are  given  in  Section  6. 


2.2.  ^oUition  of  the  Least-Squares  Problem 

After  A has  been  reduced  to  upper  triangular  form  by  the 
process  described  in  Section  2.1,  the  matrix  R in  (3)  will  be  non- 
singular (because  rank  (A)  = n)  . Tiie  unique  solution  of  the  full- 
rank  least-squares  problem  is  then  computed  as  follows.  The  matrix 
(J  (the  product  of  the  Householder  matrices  •••  H^)  reduces  A 

to  upper  triangular  form  from  the  left  and  is  orthogonal.  Because 
Euclidean  length  is  preserved  by  ortliogonal  transformation,  the 
least-squares  residual  of  the  problem  transformed  by  Q is  equal  in 
lengtti  to  the  residual  of  the  original  problem;  thus: 


IIAx  - bll,  = 


llQ(Ax  - b)!!.,  = IIRx  - Qbll  = 


wliere  tiie  vector  Qb  has  been  partitioned  into  its  first  n com- 


ponents (li)  and  rem.iinlng  (m  - n)  components  (b)  . 
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liocause  rows  n + 1 through  m of  R aro  zero,  it  is  clear 
from  inspection  of  (4)  that  tlie  residual  vector  of  the  transformed 
problem  will  be  minimized  when  the  first  n components  of  Rx  are 
equal  to  the  first  n components  of  Qb.  i.e.,  a mlnlmun!  residual 
will  be  attt'iined  for  the  vector  x that  satisfies  the  n by  n 
non-singular  linear  system: 


Rx  = b 

Tile  minimum  residual  will  have  length  Ilbll2,  since  tlie  residual  of 
the  transformed  problem  (4c)  is  zero  in  tlie  first  n components, 
by  construction  of  x to  satisfy  (5)  . 

Tlic  orthogonality  of  the  transformations  used  to  reduce  A 
to  triangular  form  is  crucial  in  solving  tlie  least-squares  problem; 
the  solution  of  the  transformed  problem  in  (4)  is  equivalent  to  the 
solution  of  tlie  original  problem  only  because  Euclidean  length  is 
preserved . 

2.3.  Summary 

In  the  full-rank  case,  tlie  linear-least  squares  problem  is 
solved  as  follows; 

(a)  ■construct  a sequence  of  Householder  transform.itions  such  that 

R 

0_ 

where  R Is  upper-triangular;  the  result  is  sometimes  described 
as  the  "QR  factorization"  of  A. 


H • • • H.H, A = QA  = 
n 2 1 
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(b)  tr.uisforni  the  rlght-hind  side,  l.e.,  form 


[ 


} n 

} m - n 

(c)  solve  the  triangular  system  Rx  = b; 

(d)  if  desired,  compute  the  residual  vector  for  the  original  problem 
by  back- transforming  the  transformed  residual: 

Ax  - b = Q^[QAx  - Qb]  = 

Step  (a)  needs  to  be  carried  out  only  once  for  a given  matrix 
A,  after  which  the.  least-squares  problem  can  be  solved  for  different 
right-hand  sides  by  repeating  steps  (h'>  through  (d)  . 


0 

b 


H • • • H.lLb  Qb  = 
n 2 1 ' 


i 

I 
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3.  11h‘  R;ink;^r)yl  i(- IcT^t  Linear  l.,east-S(7uares  Problem 

3.1.  ^on^thi  i^iu  tU'.ss  ; tjie  Mi  ninuirrM.en^t  h So  I u t i on 

If  the  matrix  A has  deficient  column  rank,  the  linear  least- 
squares  problem  is  more  complicated.  If  rank  (A)  = r,  where  r < n, 
we  can  assume  for  purposes  of  exposition  that  the  first  r columns 
of  A are  linearly  independent.  After  applying  r Householder 
transtormations  to  A on  tlie  left,  as  detcrihed  in  Section  2.1, 
tlie  resu  It  will  he  : 

n 

. Mir 

H 11,  A OA  = K = = I (6a) 

1_0J  } m - r 

where  R is  upper  trapezoidal  , since  columns  (r  ■(  I)  through  n 
of  A are  by  assc.mption  linearly  dependent  on  the  iirst  r columns 
(the  block  of  zeros  below  ii  is  absent  if  m = r) . Let  the  trans- 
formed Vv'ctor  Qh  be  similarly  partitioned: 


From  (6),  It  can  lie  seen,  as  in  the  full-rank  case,  that  a 
min  i mum- 1 ongt  li  residual  HAx  - bit.,  will  he  attained  for  any  vector 
X th.it  saMsfles  the  set  of  r equations: 


1 
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Rx  = b 


(7) 


Ik'caii.st'  R has  fewer  rows  than  columns,  the  solution  to  (7)  Is 
not  unique,  and  thus  the  solution  to  tlie  least-squares  problem  is 
not  unique. 

AmonR  all  vectors  satisfying  (7),  it  is  often  desired  to  find 
the  oni'  with  minimum  Euclidean  length,  lliere  arc  various  ways  to 
characterize  the  minimum-length  least-squares  so’utlon;  a full  discus- 
sion is  given  in  Peters  and  Wilkinson  (1970).  The  properties  of 
interest  here  derive  from  the  observation  that  tlie  r rows  of  the 
upper  trapezoidal  matrix  R are  linearly  independent,  and  span  a 
subspace  of  dimension  r in  n-space,  of  all  vectors  that  are  linear 
combinations  of  tlie  rows  of  R.  The  minimum-length  vector  that  solves 
the  least-squares  problem  is  tlie  unique  vector  satisfying  (7)  that 
lies  entirely  in  this  subspace. 

To  see  why  this  is  true,  assume  that  the  r columns  of  a 
matrix  V form  a basis  for  this  suhspace.  The  particular  representa- 
tion of  V is  not  relevant  initially  (for  example,  the  rows  of  R 
clearly  comprise  such  a basis);  additional  discussion  is  given  in 
Section  3.2.  By  definition  of  a basis,  any  vector  x in  this  sub- 
space may  be  written  as  a linear  combination  of  the  columns  of  V,  i.e.. 


X = Vw 


(8a) 
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L 


I 

I 


1 

I 


1 

T 


for  some  r-vector  w.  To  show  thot  such  an  x will  satisfy  (7), 
we  substitute  the  expression  for  x into  (7); 


RVw  = b 


(8b) 


Tlie  matrix  RV  is  r by  r,  and  is  non-singular  by  definition  of 
V as  a basis;  lienee,  w is  unique.  Thus,  when  w is  the  solution 
of  (8b),  X - Vw  satisfies  (7),  and  yields  a mlnimum-lengtli  residual 
for  the  least-squares  problem. 

To  verify  that  x has  minimum  Euclidean  length  among  all 
vectors  satisfying  (7),  we  consider  a complementary  subspace  to 
that  mentioned  above  — namely,  the  subspace  of  dimension  (n  - r) 
of  vectors  orthogonal  to  the  rows  of  R.  Let  the  (n  - r)  columns 

of  Che  matrix  V form  a basis  for  this  complementary  subspace,  so  t 

_ I 

that  every  vector  z for  wlilch  Rz  = 0 may  be  written  as;  | 

I 

1 


z = Vy 


for  some  (n  - r)-vector  y,  and 


V^\)  = 0 


Let  X be  any  vector  that  yields  a minimum-length  residual; 
then  X must  satisfy  (7): 

Rx  = b 


(<)) 
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Because  x also  satisfies  (7),  it  follows  that 

R(x  - x)  = 0 . 

Thus,  the  vector  (x  - x)  is  orthogonal  to  tlie  rows  of  R,  and,  as 
iniiicated  above,  may  be  written  as: 

X - x = Vy  , 
or 

X = X + Vy  . (10) 

Consequently,  every  vector  x satisfying  (7)  may  be  written 
in  the  form  (10).  Consider  the  Euclidean  length  of  such  a vector: 


IxH^  = II X + Vyll^  = 2^'^Vy  + HVyll^ 


_X~ 

It  follows  from  (9)  that  the  term  x Vy  vanishes,  since  x = Vw; 
^ 2 

the  expression  for  Ilxll2  then  becomes: 


so  that 


Ixll^  = ll^ll^  + IIVyl|2  , 


Ixll^  > li;i|2  . 


(11a) 


(11b) 
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Since  tlie  columns  of  V are  linearly  independent  by  construc- 
tion, equality  holds  In  (ilb)  only  when  y = 0.  Therefore,  the 
vector  X defined  by  (8a)  and  (8b)  is  the  unique  least-squares  solution 
of  minimum  Euclidean  length. 

3.2.  Reduction  of  an  Upper  Trapezoidal  Matrix  to  Upper  Triangular 
Form  by  Transformation  From  the  Right 
As  shown  in  Section  3.1,  the  minimum-length  least-squares 
solution  can  be  computed  using  a set  of  basis  vectors  for  the  subspace 
spanned  bv  the  rows  of  the  upper  trapezoidal  matrix  R.  Through 
application  of  another  special  sequence  of  Houseliolder  transformations, 
it  is  possible  to  construct  an  orthogonal  basis  for  this  subsnace 
by  rcdui'ing  R from  tji^c  r ight  to  an  upper  triangle  followed  by  a 
block  of  zeros. 

Suppose  that  there  exists  an  n by  n orthogonal  matrix  V 
that  reduces  R to  this  form,  l.e.; 

r n-r 

RV  = [R  j 0 ] , (12) 

where  R is  an  r by  r non-singular  upper  triangle. 

If  the  columns  of  V are  partitioned  accordingly: 

r n^ 

V = [V  I V ] , 


I 
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then  It  holds  that 


Hence,  the  r columns  of  V form  an  orthogonal  basis  for  the  space 
spanned  by  the  rows  of  R,  and  the  (n  - r)  columns  of  V form  an 
orthogonal  basis  for  the  subspace  of  vectors  orthogonal  to  the  rows 
of  R. 

From  the  results  in  Section  3.1,  the  minimum-length  least-squares 
solution  X may  be  written  as  x = Vw . Equation  (7),  which  specifies 
the  minimum- residual  property,  then  becomes: 

Rx  = RVw  = Rw  = b , (13) 

so  that  w is  the  solution  of  a linear  system  involving  the  triangular 
matrix  R . 

The  orthogonality  of  the  chosen  representation  of  V guarantees 

that  the  condition  number  of  the  matrix  RV  is  no  greater  than  that 

of  R;  tlius,  the  "natural"  conditioning  of  the  problem  is  not  altered 

by  an  orthogonal  V.  This  favorable  result  would  not  hold  for  some 

- -T 

other  choices  of  V — for  example,  if  V were  taken  as  R , the 

— T 

matrix  in  (13)  would  be  RR  , whose  condition  number  is  the  square  of 
cond (R) . 
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The  process  of  constructing  an  orthogonal  matrix  V that  satis- 
fies (12)  involves  the  definition  of  r Householder  transformations, 
and  is  best  illustrated  by  an  example.  Let  n = 5,  r = 3,  so  that 
R is  given  by: 


'■11 

’^12 

’'13 

*^14 

"15 

R = 

0 

*^22 

^^23 

’'24 

’’25 

0 

0 

r33 

*^34 

*^35 

The  first  step  of  the  reduction  Involves  constructing  a 
Householder  transformation  — — that  annihilates  elements  4 and 

5 in  row  3.  Since  elements  1 and  2 of  row  3 are  already  zero,  the 
corresponding  Householder  vector  will  be  non-zero  only  in  positions 
3,  4,  and  5,  and  so  will  not  alter  columns  1 and  2 of  R.  The  result 
of  applying  H,^  is  the  following: 


’’ll 

’^12 

"13 

1 

1— ■ 

‘^15 

RH^  = 

0 

'^22 

’*23 

’'24 

■^25 

0 

0 

*^33 

0 

0 

where  the  barred  elements  have  been  altered,  and 

|-  I ^ 2 2 ^ 2 ,1/2 

''■33'  = ^"33  •■34  ’■35^  • 

•r 

1 

I 
I 
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Next,  a transformation  H2  is  required  to  annihilate  elements 
•4  and  5 of  row  2.  In  order  for  the  form  of  row  3 to  be  retained,  the 
inner  product  of  row  3 and  the  vector  corresponding  to  H2  must  vanish, 
which  implies  that  the  third  component  of  the  Householder  vector 
must  be  zero.  Consequently,  the  vector  corresponding  to  H2  will 
have  non-zeros  only  in  positions  2,  4,  and  5,  and  the  effect  of  H2 
is : 


s 

s 

s 

■'ll 

*'12 

■"n 

''14 

"15 

^382  » 

0 

*^22 

''23 

0 

0 

0 

0 

0 

0 

33 

where  the  doubly  barred  elements  have  been  altered,  and 

I'  I / 2 , -2  ^-2.1/2 

^22!  - (r22  + *^25^ 

Finally,  is  constructed  to  annihilate  elements  4 and  5 in 

the  first  row;  the  vector  corresponding  to  must  have  zeros  in 

positions  2 and  3,  in  order  not  to  affect  rows  2 and  3.  The 
ultimate  result  is: 


r r r 

11  12  *^13 


*^22  *^23 

° *^33 


0 0 

0 0 

0 0 


which  is  the  desired  reduced  form. 
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Ill  summary,  the  n by  n orthogonal  matrix  V that  satisfies: 


RV  = [R  j 0] 


ts  given  by  the  product  of  r Householder  transformations,  constructed 
in  a backward  sweep  over  the  rows  as  described,  and  may  be  written 


V = H H ,•  • • H,  , 
r r-1  1 


It  should  be  emphasized  tiiat  these  transformations  do  not  have  tlie 
same  structure  as  those  used  in  the  reduction  from  the  left;  in 
general,  the  vector  corresponding  to  (applied  on  the  left)  has 

zeros  In  components  1 through  (i  - 1),  and  non-zeros  in  components 
i through  m;  whereas  the  vector  corresponding  to  (applied  on 

the  right)  has  non-zeros  in  component  1 and  components  r + 1 
tiiroug'a  n,  and  zeros  in  all  other  components. 


1.3.  ^o lution  of  the  Least-Squa re s Pro blem 

If  it  is  assumed  that  the  first  r columns  of  the  matil.<  A 


are  linearly  independent,  the  minimum-length  least-squares  solution 
may  be  computed  as  follows: 
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(a)  Reduce  A to  upper  trapezoidal  form  by  application  of  r 


Householder  transformations  on  tlie  left,  using  the  procedure 
described  in  Section  2.1,  so  that 


H 

r 


H^A  £ QA  = 


I (b)  Transform  the  right-hand  side  by  applying  Q,  i.e.,  form 


•i 


Qb 


b 


} r 
} m-r 


(c)  Construct  a sequence  of  r Householder  transformations  to  be 

applied  on  the  right  to  reduce  R to  upper  triangular  form,  as 
described  in  Section  3.2,  so  that 


r n-r 

QAH^  •••  H^  QAV  = 

where  V is  the  product  of  the  second  set  of  Householder  trans- 
fo  niia  Lions. 
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(d)  Solve  the  non-singular  upper  triangular  system 

Rw  = b 

(e)  Transform  w by  applying  V,  the  first  r columns  of  V, 
yielding  the  solution  x as 

X = Vv 

If  tlie  minimum-length  solution  is  not  required,  a least-squares 
solution  ran  be  computed  by  carrying  out  steps  (a)  and  (b),  followed 
by: 

(c')  partition  R as  follows: 

r n-r 
[R  S], 

so  that  R is  the  left-most  r by  r submatrix;  R must  be  non- 
singular if  the  first  r columns  arc  linearly  Independent. 

(d')  solve  kx  = b,  and  let 
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Tlie  solution  obtained  in  this  way  will  satisfy  equation  (7), 


ri 


I 


and  is  the  least-squares  solution  of  the  problem  corresponding  to 
the  first  r columns  of  A and  the  given  right-hand  side.  In 
general,  it  is  not  the  minimum-length  solution  corresponding  to  all 
n columns  of  A. 
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4.  Estimation  of  Rank 


r 

f 

[ 


I 

t 

P 


4.1.  J^i_t_f  ten  l^^i  t^s 

The  procedures  described  in  Section  3 for  the  rank-deficient 
least-squares  problem  assumed  that  the  rank  of  A was  known  a priori 
to  bo  r,  and  that  the  first  r columns  of  A were  linearlv  indepen- 
dent. Tliese  assumptions  are  not  realistic,  and  were  Introduced  only 
for  simplicity  of  presentation.  In  attempting  to  solve  a least- 
squares  problem  whore  the  matri.x  A is  of  unknown  rank,  it  is  obviously 
necessary  to  estimate  the  rank  during  the  course  of  the  computation, 
and  thereby  to  determine  a set  of  linearly  independent  columns. 

Unfortunately,  the  definition  of  "rank"  in  the  context  of 
computation  with  floating  point  arithmetic  and  inaccurate  data  depends 
on  the  problem.  The  question  can  never  be  resolved  in  a specific 
case  without  making  an  explicit  judgment  about  the  scaling,  i.e.,  a 
determination  ns  to  which  quantities  can  be  considered  as  "negligible" 
(for  an  excellent  discussion,  see  Golub,  Klema,  and  Stewart,  1976). 

As  an  illustration  of  the  complexity  of  the  issue,  consider 
the  matrix 

'1  r 

_0  c_ 

where  £ Is  not  zero,  but  is  small  relative  to  unity.  Mathematically, 
the  two  columns  are  linearly  Independent,  and  the  matrix  has  rank  2.  . 


! 


( 


i 

! 


I 
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In  practice,  however,  the  second  vector  may  be  a computed  version  of 
the  first,  so  that  numerically  the  two  columns  should  be  considered 
"equivalent",  even  though  one  is  not  an  exact  multiple  of  the  other; 
in  this  event,  the  matrix  has  "rank"  1.  Thus,  the  decision  as  to 
whether  these  two  vectors  are  linearly  independent  depends  on  whether 
or  not  the  value  of  c is  "negligible";  even  in  the  simplest  2 by 
2 case,  the  issue  of  rank  can  be  resolved  only  by  considering  the 
nature  of  the  underlying  problem,  and  the  origins  of  the  data.  For 
larger  values  of  n,  a satisfactory  automatic  technique  for  determin- 
ing the  rank  of  a matrix  is  still  more  complicated. 

With  exact  arithmetic,  linear  dependence  among  a set  of  columns 
would  reveal  Itself  during  the  reduction  to  upper  triangular  form 
described  in  Section  2.1.  If  the  (k  + l)-st  column  were  a linear 
combination  of  the  previous  k columns,  components  k + 1 through 
m of  the  transformed  dependent  column  (the  "remaining  column")  would 
be  exactly  zero.  With  finite  precision  computation,  one  might  accord- 
ingly hope  that  tlie  norm  of  a remaining  column  would  be  "small"  if 
that  column  were  "nearly"  linearly  dependent  on  the  previous  columns. 

Unfortunately,  this  hope  is  not  realized,  since  it  is  equivalent 
to  expecting  that  an  ill-conditioned  triangular  matrix  will  have  at 
least  one  small  diagonal  element.  Triangular  matrices  exist  that 
have  no  "small"  diagonal  elements,  yet  are  arbitrarily  badly  conditioned  — 
for  example,  the  famous  n by  n matrix: 
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I 
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whose  rondition  number  is  of  order  d”  and  which,  even  for  moderate 
values  ol  n,  could  reasonably  he  termed  "numerically  rank-deficient". 

Tlie  Implication  of  this  example  for  the  estimation  of  rank 
during  triangular  reduction  is  clear:  if  one  were  reducing  to 

upper  triangular  form,  no  transformations  would  be  necessary  (since 
it  is  already  an  upper  triangle),  and  no  indication  of  t lie  near  linear 
dependence  among  the  columns  would  be  given.  Consequently,  a strategy 
based  on  the  expect.ation  of  a "small"  remaining  colun:n  in  th-  presence 
of  near  linear  dependence  can  not  be  guaranteed. 

4.2.  Column  Interchange  Strategy 

Despite  tlie  rather  artificial  example  given  in  Section  4.1, 
numerical  linear  dependence  is  almost  Invariably  revealed  in  practice 
by  a "small"  remaining  column  during  the  reduction  from  the  left  (see 
further  remarks  that  strengthen  this  observation  in  Golub,  Klema,  and 
Stewart,  1976).  Thus,  an  obvious  strategy  for  estimating  rank  and 
selecting  a set  of  linearly  independent  set  of  columns  is  to  carry  out 
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column  intorcliangos  clurlnj;  tlie  reduction  from  the  left,  and  to  reduce 
at  tlie  k-th  step  the  remaininn  column  of  "largest"  norm,  in  order  to 
move  tlie  "small"  columns  to  the  end.  This  strategy  is  discussed  in 
Peters  and  Wilkinson  (1970)  and  Lawson  and  Hanson  (1974),  and  is 
sometimes  called  "column  pivoting".  It  is  ni^t  guaranteed  to  move  tlie 
best-conditioned  set  of  columns  to  the  "front"  (left)  of  the  matrix, 
but  in  general  it  tends  to  allow  the  "mo.st  Independent"  columns  to  bo 
processed  first.  The  effect  of  interchanging  columns  is  to  Introduce 
a permutation  matrix  on  the  right  of  the  matrix  A in  (6a);  the  hoped- 
for  configurat ion  then  becomes: 


H • • • H , AP  QAP  = R 
r 1 


(14) 


where 


P is  a permutation  matrix. 

In  practice,  the  computed  version  of  (14)  will  be 


n-r 


QAP 


m - r 


where  llpll  is  "small".  A detailed  discussion  of  the  significance 
of  L is  given  in  Colub,  Klema,  and  Stewart  (1976). 
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In  solving  tho  least-scjuares  problem  with  a column  intercliange 
strategy,  steps  (e)  and  (d')  of  the  procedures  given  in  Section  3.3 
must  be  expanded  to  include  application  of  tiie  permutation  matrix  P 
to  the  solution  x (i.e.,  appropriate  components  of  x are  interchanged) 

4.3.  Oonsiderat ions  in  Estimating  the  Rank 
4.3.1.  "Size"  of  a column 

Tile  idea  of  estimating  rank  witli  a column  interchange  strategy 
to  reduce  the  "largest”  remaining  column  at  each  step  leads  to  an 
interesting  question:  what  sliould  be  the  definition  of  the  "size" 
of  a column?  The  most  obvious  measure  of  the  size  of  a vector  is 
simply  its  Kuelidean  Length;  but  since  the  Houseliolder  reduction  from 
the  left  acts  on  each  column  independently,  this  definition  is  appro- 
priate only  if  one  assumes  tb.at  the  matrix  to  be  reduced  has  bei  n 
scaled  initially  so  that  the  same  definition  of  "small"  applies  uni- 
formly to  all  columns.  In  many  case's,  however,  c'ach  column  represents 
observations  of  a particular  variable,  and  the  values  of  the  variables 

may  not  display  a constant  scaling  with  respect  to  one  another  — for 

_2 

example,  the  average  magnitude  of  one  variable  might  be  10  , while 

another  might  be  of  order  10^.  In  such  an  instance,  it  would  be 
wrong  to  reg.ird  tlu'  column  correspond  i ng  to  the  first  variable  as  seven 
orders  of  m.ignitude  "smaller"  than  the  column  corresponding  to  the 
second  variable.  Clearly,  it  is  necessary  to  exercise  care  in  the 
definitio\i  of  "size"  in  older  to  use  a column  intercliange  strategy 
to  estimate  rank.  This  topic  will  be  discussed  in  more  detail  in 
Section  6.4.11. 
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A.  3. 2.  Factors  inf  1 uonclnf,  the  estimate  of  rank 

In  many  applications  it  makes  practical  sense  to  err  on  the 
side  of  ni^de res ^iana ting  the  rank,  for  two  reasons.  First,  a lower 
estimate  of  rank  often  yields  a more  reasonable  (l.e.,  smaller  in 
norm)  com])nted  solution,  but  does  not  appreciably  increase  the  size 
of  the  residual.  Peters  and  Wilkinson  (1970)  cite  the  followinR 
example: 


*6  3.0000  OOOOO' 

~3.0000  OOOOO" 

A = 

A 1.9999  99998 

h = 

2.000A  00000 

_2  1.0000  0000 3_ 

0.999A  00000 

If  the  matrix  A is  considered  to  have  rank  one,  the  solution 

T -A 

is  very  close  to  (O.A,  0.2)  , and  the  residual  vector  is  of  order  10 

On  the  other  hand,  if  tlie  matrix  is  considered  to  have  rank  two,  the 

I 5 

solution  will  he  very  close  to  (10  ,-2  x 10  +1),  and  the  residual 

-9 

is  of  order  10  . Altliough  the  very  large  solution  vector  yields 

a smaller  residual,  such  a solution  may  be  nonsensical  in  practice, 

and,  in  fact,  the  first  vector  is  almost  certain  to  be  a much  "better" 

solution,  despite  the  larger  residual. 

A second  reason  for  tending  to  favor  a conservative  estimate 

of  the  rank  is  that  the  perturbation  analysis  of  the  least-squares 

problem  shows  that  the  relative  change  in  the  exact  solution  can  include 
2 

a f.actor  (cond(A))  for  an  incompatible  right-hand  side  (see  Stewart, 
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1973).  Thus,  in  order  to  avoid  an  ill-conditioned  problem,  one  j 

miglit  wish  to  accept  a smaller  value  of  tlie  rank,  and  thereby  to 

select  a set  of  columns  that  are  presumed  to  be  "strongly"  linearly 

independent.  Although  a larger  estimate  of  the  rank  would  be  possible 

with  a more  stringent  criterion  for  dependence,  the  conditioning 

of  tlie  least-squares  problem  might  then  reflect  the  square  of  a much 

increased  condition  number  for  the  extended  set  of  columns. 

There  are,  nonetheless,  instances  where  one  may  wish  to  allow 
a very  liberal  estimate  of  the  rank  to  be  made  during  tlie  triangular 
reduction.  A notable  example  occurs  for  linearly  constrained  least- 
squares  problems,  where  the  linear  constraints  are  often  known  to 
remove  the  difficulties  with  near  linear  dependence  in  the  unconstrained 
problem.  In  this  case,  it  would  not  be  appropriate  to  terminate  tin- 
triangular  reduction  prematurely,  since  the  possible  dangers  from 
over-estimating  the  rank  would  be  avoided. 

In  summary,  the  best  strategy  for  estimating  the  rank  of  a ' 

matrix  during  reduction  to  triangular  form  depends  on  properties  of  1 

the  ultimate  problem  to  be  solved.  j 


1 

J 


f 
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I 5.  The  Complete  Orthogonal  Factorization 

! The  column  Interchange  strategy  described  In  Section  4.2 

allows  the  procedures  of  Sections  2.1  and  3.2  to  be  applied  to  compute 
[ the  complete  orthogonal  factorization  of  a general  real  m by  n 

matrix  A of  rank  r.  If  column  interchanges  are  carried  out  to 
j ensure  that  r linearly  independent  columns  are  processed  first,  the 

I 

I.  configuration  after  the  reduction  from  the  left  is; 


QAP 


where  R is  upper  trapezoidal,  and  P is  a permutation  matrix  that 
specifies  the  column  interchanges.  The  reduction  of  R to  upper 
triangular  form  is  then  carried  out  if  necessary  (see  Section  3.2), 
and  tlie  orthogonal  permutation  matrix  P is  absorbed  into  the  orthogonal 
matrix  arising  from  the  reduction  on  the  right  (this  simply  inter- 
changes the  rows  of  the  latter  matrix).  The  final  result  is 


n-r 


QAV 


R 0 
0 


} m - r 


and  the  complete  orthogonal  factorization  of  A is  given  by  the 
matrices  Q,  V,  and  R. 
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This  decomposition  is  extremely  useful  in  many  applications. 
Let  Q be  partitioned  as: 


} r 


Q2  } m - r 


'Hie  columns  of  form  an  orthogonal  basis  for  the  space 

spanned  by  the  columns  of  A,  and  the  columns  of  ©2  form  an  ortho- 
gonal basis  for  the  set  of  vectors  orthogonal  to  the  columns  of  A. 
In  addition,  the  last  (n  - r)  columns  of  V form  an  orthogonal 
basis  for  the  set  of  vectors  ortiiogonal  to  the  rows  of  A. 

Any  m-vector  y may  be  written  as  the  sum  of  two  orthogonal 
parts,  which  lie  respectively  in  the  column  space  of  A and  it.s 
orthogonal  complement,  i.e.. 


y = + Q2y2  ■ 


The  vectors  y^  and  'j ^ can  be  obtained  by  forming  Qy  since; 


(Qfyi  + Q2y2) 


Tims,  the  matrix  Q can  be  applied  to  compute  the  projection 
of  a vector  into  the  column  space  of  A or  its  orthogonal  complement, 


i 

t 


without  using  projection  matrices  (see  Stewart,  1973,  for  further 
details).  Since  Q is  tlie  product  of  a sequence  of  Householder 
transformations,  it  is  not  necessary  for  this  purpose  to  form  Q 
explicitly,  but  merely  to  apply  the  transformations  to  the  given 
vector . 

The  bases  represented  by  the  columns  of  Q are  also  useful 
in  other  contexts  — in  particular,  for  solving  optimization  problems 
with  linear  constraints.  Suppose  that  a problem  contains  the  n 
linearly  Independent  linear  constraints 

A%  = c . (15) 

wliere  y is  an  m-voctor.  The  complete  orthogonal  factorization  of 
A provides  a straightforward  way  to  reduce  the  dimensionality  of 
the  optimization  problem.  A vector  v that  satisfies  (15)  can  be 
found  by  solving  the  least-squares  problem: 

min  llA^y  - cH^  • 

which  will  be  compatible  if  A has  full  rank.  Every  vector  y 
satisfying  (15)  can  then  be  written  as: 

y “ y + Q2''  • 
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for  some  (m  - n)-vecLor  v,  because  the  columns  of  Q2 
orthogonal  basis  for  the  set  of  vectors  orthogonal  to  the  columns 
ot  A.  Therefore,  only  the  optimal  choice  of  v is  unknown,  and 
I any  further  minimization  will  occur  with  respect  to  these  (m  - n) 

variables.  Furthermore,  every  iterate  y can  be  represented  in 
the  form  (16),  and  will  satisfy  the  linear  constraints  by  construction 
• of  Q,^.  Tliese  ideas  may  also  be  applied  in  other  contexts  — for 

I example,  linear  Inequality  constraints  (Gill  and  Murray,  1974),  or 

I separable  nonlinear  least-squares  with  separable  nonlinear  equality 

j constraints  (Kaufman  and  Pereyra,  1978). 
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6.  Doouiuiyitation;  Compu^iticmal  Details 

There  are  three  major  catejtories  of  subroutines  in  this  collec- 
tion, in  iiddition  to  the  general  subroutine  MNLNLS: 

(1)  subroutines  to  carrv  out  the  Houseliolder  r“duction  of  a given 
matrix  to  suitable  form.  Tins  group  includes  the  subroutines 
HKKDL.  HMUl.TC,  I.KNSQ,  and  HRF.DR; 

(2)  subroutines  that  solve  die  least-squares  problem  after  the 

matrix  has  been  factorized.  This  group  includes  the  subroutines 

QRVSLV.  QMIT.VC,  VNiLT.VC,  TRSLV,  TRTSLV,  and  UNSCRM; 

(T)  subroutines  lliat  explicitly  form  the  orthogonal  matrices 
T 

Q(or  Q ) and  V.  This  group  includes  the  subroutines  FRMORT 
and  FORMV. 

The  subroutines  in  this  package  have  been  designed  to  include 
substantial  flexibility,  to  allow  their  use  in  widely  varying  contexts. 
The  user  who  simply  wants  to  solve  the  linear  least-squares  problem, 
and  lines  not  wish  to  be  concerned  with  any  details  of  the  algorithms 
or  fine  points  of  the  code,  should  use  the  subroutine  MNI.NhS;  MNLNLS 
is  documented  in  Section  b.7,  and  auto.iuitically  calls  the  appropriate 
routines.  Otliorwise,  tlie  following  brief  alphaiietical  guide  indicates 
the  purpose  of  the  subroutines: 

— FORMV  (Section  6.1)  — form  the  explicit  ortliogonal  matrix  V; 

— FRMORT  (Section  6.2)  — form  the  explicit  ortliogonal  matrix  Q 
(or  Q*^); 

14 


— IIMIU.'I'C;  (Scctiini  {).  ^)  — construct  .iiid  apply  a single  lloiisclin  1 dor 
t tans  format  ion  to  a set  of  columns; 

— lIRKDL  (Section  6.4)  — carry  out  the  Householder  reduction  from 
the  left,  including  estimation  of  rank; 

— HRICDR  (Section  6.5)  — carry  out  the  Householder  rt^duction  from 
t lie  rii'.ht  of  an  upper  trapezoidal  matrix; 

— LENSQ  (Section  6.6)  — compute  tlie  squared  Kuclidean  length  of 
a vector; 

— MMl.NI.S  (Si'Ction  6.7)  — compute  the  mini  mum- length  least-squares 

I 

sol ut  ion ; 

I 

— QMUI.VC  (Section  6.8)  — apply  to  a vector  the  sequence  of  transfor- 
mations constructed  t>i  reduce  a matrix  from  the  left; 

— QRVSLV  (Section  6.9)  — solve  the  least-squares  problem  after  tlie 
matrix  has  been  reduced  to  an  appropriate  form; 

— TRSl.V  (Section  ti.lO)  — solve  a triangular  system  of  linear 
equations ; 

— TRTSIA'  (Section  6.11)  — solve  the  transpose  of  a triangul.ir  system 
of  linear  equations; 

— UNSCRM  (Section  6.12)  — apply  a permutation  to  a vector; 

— V.Ml'I.VC  (Section  6.13)  — apply  to  a vector  the  sequence  of  transfor- 
mations construi'ted  to  reduce  a trapezoidal  matrix  from  the  right. 
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b . I . Siibrc^tJjio  R)R.W 
f) . 1 . 1 . Purpose 

The  subroutine  FOKMV  forms  the  explicit  n by  n orthogonal 
matrix  V that  is  the  product  of  r Householder  transformations  applied 
on  the  right  by  the  subroutine  HRKDR  to  reduce  an  r by  n upper 
trapezoidal  matrix  R to  an  r by  r upper  triangle  followed  by 
a block  of  zeros.  The  first  r columns  of  V form  an  orthogonal 
basis  for  the  set  of  vectors  spanned  by  the  rows  of  R;  the  last 
(n  - r)  columns  of  V form  an  orthogonal  basis  for  the  set  of  vectors 
orthogonal  to  the  rows  of  R.  Details  are  given  in  Section  3.2. 

6.1.2.  Description  of  method 

The  matrix  V is  the  prodvict  of  tlie  r transformations  applied 
to  R on  the  riglH:; 


where  the  vector  corresponding  to  H,  has  non-zeros  in  components 
j,  and  (r  + 1)  titrough  n.  To  construct  V,  those  transformations 
are  multiplii'd  together  beginning  at  the  right,  after  is  formed 

explicitly.  Finally,  any  column  interclianges  carried  out  during  the 
reduction  from  the  left  to  upper  trapezoidal  form  are  incorporated 
by  applying  tlie  appropriate  permutation  matrix  on  the  left  (inter- 
ch.mging  the  rows  of  the  orthogonal  matrix). 
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Keywords 

Householder  reducCiou  from  right;  complete  orthogonal  factor! 

zation . 

6.1.4.  Source  language 

Fortran.  Tiie  code  in  FORMV  has  been  checked  by  the  PFORT 
verifier,  and  is  WAi'FIV-compatiblc.  All  variables  and  functions  are 
explicitly  declared. 

6.1.5.  Specification  and  parameters 
See  accompanying  listing. 

6.1.6.  Error  indicators 

See  accompanying  listing  (the  description  of  the  parameti 
LKRROR) . 

6.1.7.  Auxiliary  routines 

FORMV  calls  the  standard  function  DABS. 

6.1.8.  Program  size 

71  Fortran  source  statements. 

6.1.9.  Array  storage 


No  locally  declared  arrays. 


6.1.10.  Timing 


The  number  of  arithmetic  operations  required  to  form  V Is 
2 2 

ot  approximate  order  r (n  - r)  + 2r(n  - r)  . 


6.1.11.  Further  Comments 


I 


I 


StJHFCUriNs:  F0?“'/  I NCJL,  NDI*',  gpV,  HVFCa,  IPFffl. 

1 N/DIM,  'c‘T,  irBFOa  ) 

c 

:n''EgE(»  nc-l,  NPrNK,  NVDir,  lebfor 

rN"’EGE:<  !PF“'1(N.'ni 

DC'rP’.i:  t'BECISION  Ol-  V C."^  Ih  , 'JCw'D,  HV7v~  F (SCOL)  , V(NVDr!1,  NCO  LJ  , 

1 fSE^L) 

r 

c 



c 

c 

C THF  jlJbSOU'’  i FCP  .V  IE  Uf»FD  IN  CO  NJ  UN  C'^T'^N  WITH  THE  SUB  R OU ’’'I  NSS 

C HPIDL  AND  HEtOf-  ID  Ci-.Ft'IT:  fSRT  OF  ThF  CIKFLETF  CSIKOOCNRL 
C FACT0FI7A::CN  .F  a CrN'^rAL  FEAL  “ATFiy. 

C 

C FOr-  EVERY  !■  :;a  I NE'W  EY  NCOT.  SEAI  HATPIX  ORV  OP  PANE  VPANK, 

C '"tlPFE  FXIjT  an  tiH',  * :Y  NFOW  "FEHOr,  )NAL  ."'.TFrY  C,  AN  NFU  PY  NCOL 
C O^THOOrjfO  FAIFiX  V,  ANC  AN  NRAKK  3Y  NPANK  U F F t F I F I A N'l  UL  AP  •'.ATNIV 
C P,  SUCF  riiA.  * 
r 

C C*grV*V=<P  o> 

c < 0 > . 

c 

r "HP  1ATHIX  0 IS  '^HO  FPCDUCT  OF  NPANK  HruSEHCLPEP  OR  A NS  F OK  .“A  T I iN  S 

C CONSTRUCTED  BY  IHc  SLPT  'IlINF  HPEPL.  TF  NFANK  . ” C . NfOI,  THF  HA.TUX 
C V is  FT'':’Ly  IHE  rFFrUTAIICN  “AT?  IX  DbFINFD  lY  ANY  COLUF.N  INTFF- 
C CMANC-FS  X.AD;:,  OUPINO  rXECOriOK  ^F  HFEDL.  IF  NPANK  .LT.  NCJL,  V IF 
C THE  F'-''I»nC?  IF  l.iE  ? I F " t' 1 AI  1 ON  YATPIX  A FD  A S'’Cti'^NCS  CF  NPANK 
C H'USEFC  LjER  TRANSFCFF  a::  'NS  CCNSirUCIEO  DUFINS  FXFC'IUjN  OF  THF 
C S’lBFOUTINE  HPXDI  I.AFT  ( N CO  L -SR  A N K ) CClijyNS  CF  V POPP.  AN 

C Ci^rHOGCNA^  -MjIS  F F '.Hr  SUTSPaCF  0^  VPCTOFS  0''TH''G'>NA  L TO  TH'’ 

C FCWS  'F  THE  -.iiIGINAC  HATRIX  OFV. 

C 

c 

c CHr  JEO'IENI  F '■?  K U;.  E H.;  LOE  P TF  A NS'i' CP"  A I I ^ N E 'PPLIFr  ON  '"H'’ 

C - lOHT  PY  Hi  k..OP  .v\v  3;/  W'-'ITIFN  AS 

C ? (NPANK)  * P (NPA  N*'"- 1 ) ♦ . . . ♦ Pf2)  * ?{U, 

C 

C WHFPF  F(J)  TO  A HIU.  FHCLDEP  KATFIX  CKCSEN  T ' FEDUCF.  PC  >.  0 

C THE  AFPRoPtlAiS  r.F<.  TUP  VECTCP  U(J)  O?  F F .SPC  N 01  NO  TO  P(0)  IS  0=“ 

C THE  FCLIOwING  SP'^C'aL  FOP.":  COFFONENT''  1 IHFO'irF  fJ-1)  APE  FFP^. 

C CC"FONEVr  J IS  NvN-Z.KC,  CC".rON''NTS  (J»  1)  ''1RCU0H  NPANK  A"F.  ZERO, 

C AND  CCHFO  Jr  NT.,  N.RANK*!  r(''''r'nH  NCOL  APT  NON-Z"PC.  rUFING  EXECUFION 
C OF  FOFFV,  'iflFS'’’  * F NS  FO?  NA  I ^ C NS  AFF  HULTIFLI-r  '"''O^'TH'^F  FPOF  PISH'" 

C TC  LEFT,  .AKING  ADVANTAGE,  ' K THE  SPECIAL  STFUCLUFE  OF  HU'  U{J)  TO 
C SAVF  A.PITH.X^riC  jlT.  R A"!  0 . THEN,  THF  PFF"UTATT''N  HATFIX  THAT 
C DFFINtS  ANY  CvLU-IN  I.'.  I E^’CHA  NG^S  HADE  D'lPING  THE  F EDUCTION  FPJ" 

C LEF'”  IS  APPLIED  TO  YIELF  '’’HF  DESIK’FD  MATRIX,  V. 


1 

I 
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n I'l  r,  o o r.  o o c o n o n 


i 


I 


C roPMV  CALLS  1 H.P  SURHOUTINE  'INSC^fl. 

C 

c 

C THU  FOR''AL  PAEAFETEFS  FOFF-V  APE  : 


C 

V. 

t 

c 

c 

V. 

c 

r 

c 

v« 

C 


NCOI  -- 

TNTS3EK,  INPUT  ONLY. 

THE  NU.-BEP  wF  CCLUKNS  CF  THE  OPIRINAL  (1ATRIX,  QPV. 


NDT!*  -- 

INTE.iER,  INPUT  ONLY. 

TH'’  DLCIAKED  POW  DIMENSION  OF  THE  MATRIX  QRV  IN  THE  CALLING 
PUH-PFU  GPAM. 

OPV  — 

roiJULL  PRECIS  ILN  APFAY,  OF  PECLAREC  DIM3NSICK  NPIM  BY  NCCL, 
!NPU.'  ONLY. 

THE  lAxnlX  Jl-  V C'JPR  ESPCNDS  TO  THE  CPICINAL  MATRIX  THAT  WAS 
^EDUCED  OY  HREDL  AND  HFFDR.  ITS  CCN'^EMS  MUST  NCI  DE  ALTERED 
AF’^ER  EXIT  PR  CM  HFEDR  AND  BEFORE  ENTRY  TO  FCPMV. 


NEANK  -- 

INTEGER,  INPUT  ONLY. 

;HF  NUMERICAL  RANK  '^F  THE  OFIRINAL  MATRIX  QRV,  AS  DFTERMINED  BY 
HPEDL.  THE  VALUE  CF  NFANK  MUST  NOT  BF  ALTFFFD  AFTER  EXIT  FPOM 
IIFEDL  AND  EEFJRF  FN  IP  Y TO  FOR  MV. 


HVECP  -- 

COUBi.t  PPECIsrON  VECTOR,  OF  LENGTH  NCCL,  INPUT  ONLY. 

THE  VECTOR  MVSCF  IS  GENEPATEO  BY  HREDF,  AND  ITS  CONTENTS  MUST 
C NOT  JE  ALTERED  AEIF.P  EXIT  FROM  HPEDR  AND  REF''FE  ENTFY  TO  FORMV. 

C 

C IPEPM  -- 

C INTEGER  VEl.OR,  F LENGTH  VCd,  INPUT  CNLY. 

C THE  VFCTCF  IPERM  IS  G F.N  HP  A'"  ED  BY  HPEDI,  AND  ITS  CONTENTS  MUST 

C NOT  BE  ALTERED  APIF’^  EXIT  FROM  HREtL  AND  BEFCPE  ENTRY  TC  FOPMV. 

C 

C NVDIM  -- 

C INTEGER,  INPUT  ONLY. 

C THE  UECIAREO  RCW  DIMENSION  OF  THE  MATRIX  V IN  THE  CALLING 

C EUB-PR^ GRAM.  MUST  BE  .GE.  NCOL. 

C 

C V - 

c 
c 
c 
c 
c 

C TEMP  -- 
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rCUtiLE  PRECISIcS  AFPAY,  OF  CONCEPTUAL  DIMENSICN  NCPI.  BY  NCOL, 
ANT  DECLARED  DIM'’NSICN  NVDIM  BY  NCOL.  OUTPUT  ONLY.  ON  EXIT 
FPOM  FORMV,  THE  V *’ A . F 1 X ONTAINS  THE  DESIRED  ORTHOGONAL 
■'AThlX,  V. 


DOliaLF  EFFCrSICN  V’^'C^'CF,  1F  LENGTH  NCCl,  CHTEDT  ONLY. 
USED  FOR  IN.'ES'ZI::,^  :F  SlJRAGE  DURING  FXECUTICN  OF  FOR.*V, 

LEPFOF  -- 

INTEGER,  OUTPUT  ONLY. 

AN  ERROR  INOICAT  P,  WITH  THE  FJLLTWING  P^S.STFLE  VALUES: 
LSRRJh  = 0 : NO  ERROFf,  NOEFAL  T S F 11  N A TIO  N . 

LEFRJE  = 1 : INVLID  INPUT  A R A Ft’ T '=' R . 


’*•*  AUTHJF3: 


.-ARGARET  H.  WEIGHT,  STEVEN  C.  GLAES-AN 
SYSTEM  OPTir  IZAT-ION  LABORATORY 
DiPART:;FNT  OF  OPERATICNS  RFSEAPCH 
Si’ANFOPD  UNIVERSITY 
SIANE’PR,  CALIFORNIA  941C5 


pa: 


DSCKHBFP  1977 


DECIAPATION  CP  LOCAL  VAFJAPLFS. 

INTEGER  I,  J,  JE1,  K,  NCCL"!,  NPNKPI 
DOUBLE  PRECISION  BETA,  DOT,  PRDO,  EROD 

CECLAFATICN  OP  STANDARD  FUNCTIONS. 

DOUBLE  PRECISION  DABS 


TEST  FCP  INVALID  VALUES  -F  NPAVK  AND  NCCL. 


IEEP'’^F  - 1 

■'F  (NCOL  . LE  . 0 .01.  NPANK  .LF.  0 .OR.  NRANN  .GT.  KCOL)  PETrjRN 
DO  20  J ^ 1,  NCCL 

DC  10  I = 1,  NCOL 

V (I, J)  = 0, 0C*0 
10  CONTINUE 

V (J,  J)  = 1 . or*c 
2 0 CONTiN’Jt 
LFIF^R  = 0 

IF  (NRANK  .FC.  NC'*  L)  GO  TO  2 10 
NPNKPI  = NRANK  ♦ 1 


OoOTnorjf'Jri'^o  cii".  ri  r-^oor; 


C 

c 

FCPK  THE  alGHT-HT.?!  HOU,S  iHCI.DE  R rATRIX  fXPLICITLY. 


TF  ((’VECt«(1)  .£0.  O.OD  + 0^  TO  TO  70 

BETA  = 1 .CD  + O/C-ABS  {bV2Cr<  ( 1)  ) 

PR'^D  = l.OD+C 

IF  (HV£CF(1)  .LT.  O.OD  + 0)  ORGD  = -1.00*0 

V(1,1)  = 1.00*0  - PF  ^0*HVECR  ( 1) 

DC  3"  J = -JESKI’  1,  NC-'I. 

V f 1, J)  = -EF0D*0FV  (1 ,J) 

V(J,  1)  = V (1,0) 

30  CONTI.S.'IE 

IE  (NRNKPI  .FO.  N^CLI  Oe  TO  ftO 
NCCL.r  1 = NCOl  - 1 

DO  50  J = NPNKPl,  NCOLJII 
F'^DJ  = V ( 1 ,J) 

V (j,j)  = 1.C0*''  - ':‘Poo*Gn  V ( 1 ,j) 

JP1  = J ♦ 1 

DC  4 0 I = Jr1,  *!C^L 

V (1 ,0)  = -rpDj*oP  V ( 1 , 1) 

V ( J,  1)  = '/  (I,  J) 

CON  TIN  HE 

5C  CONTINUE 

FILL  L.N  TFE  (NcjL,  NJvLI  FL^MEVT. 

5C  CONTI  SUE 

V(NCCL,NCOL)  = 1.CL*0  - FF  T A*0F  V ( 1 , N CO  LI  « * 2 
70  IP  (N’HANK  .EO.  1)  R'jrOPN 


I 

i 


•X  MI  TILL  y TOGHi’HfP  IHF  b^’rAIVINT  (XHANK-I)  I R AN  S FC  P M AT  IONS  , 
APPIVENG  EACH  NEW  TP  A N S F.)P  M.  A T TO  N FPCF  THE  LFFT. 


DC  2 00  K = 2,  NPA^'K 

PPTA  i;  THE  NTFFALT7.It:G  FACTOP  FOP  TH^  F-TH  TPA  NSFC  P"!  ATI  ON 

IF  (HVECP(K-)  .FO.  0.0r*0)  r.o  TC  200 
PFTA  = 1.C.)*0/OAtS  (HVrCF  (K)  1 
C 

c C'T.Y  ROW  K and  R ^ WS  NPANK*1  THPOnCH  NCOL  OP  IH  F CHPRENT  V 

c'  APL  aliepeo  by  application  cf  thc  k-th  tfafsfofh ation  on 

C LEFT. 


I 


C COLUMN  K C?  TH'-  AI.’.EtVO  V SIMPLE  PFrOH^IS  TPF  K -TH  COLU.*!N  OF 

C Ti'i  K-rH  HOIJ  JSriCLOEB  T BANS  FOB"  ATI  ON  . 

C 

PPOD  = l.CTtC 

IF  (HVECF(.K)  .IT.  0.0r+0»  PROD  = -1.0n+0 
V(K,K)  = A. 01-0  - PROD’HVrcP  (K\ 

D3  ICC  r = NPWKPI,  NCCL 
7(1, K)  = -PFCD*gBV  (K,IJ 
100  CONriN'JE 
C 
C 

C ALTEF  CCMIINS  1 IHPOOnH  K-1  AND  NPANF  ♦ 1 THFOtTGH  NCOL, 

C USING  A SlNGlf  D-  i-CCr  FO"  COMPACTNESS. 

C 

C THF.  LOTP  T SIAXc'r.  CNT  130  IS  A LCCP  CVEP  THU  COLUMNS  OF  THE 

C CUHNf-M  V. 

C 

DC  130  J = 1,  hCCl 
C 

C OG  NOTHIN':  IF  J .or.  K AND  J .IS.  NSANK.  THESE 

C JCLUMNS  ape  I'NALTEFFD  EFCAltSE  THFIP  TNNEP  FECDUCT  WITH 

C THL  K-TH  H;u;oWCLDEE  VECTOR  IS  /^SPO. 

C 

IF  (J  . ;E  . K .AND.  J .LE.  NPANA)  GC  ir  130 
C 

c 

C FCB';  THE  INN’"’  PRCnUCT  '’F  THE  K-TH  HCU'EHCLDEP  VECTOR 

C iND  IH?  J-i;i  COLUMN  OF  THE  CURFFNT  V. 

r- 

DOT  = 0.0  0*0 

01.  110  1=  NRNKPI,  NCOl 

DCT  = DOT  * UR  V ( K , I)  ‘V  (I , J) 

1)0  JCNi'lNU.E 

??.0  = OFTA^DOT 
C 

C ALTER  FIU  K tY  SDPTP  AC'’' ING  A MllITTFLE  CF  THE  K-TH  ELEMENT  OF 

C THE  HCU:LHCLDER  VECTOR  ( STCRED  IN  HVECP(K)  ). 

C 

V(K,J)  - V(S,J)  - FFJO’HVECR  (K) 

C 

C ALTER  ROWS  N^ANK*)  THROCGH  NCOL  OF  COLUMN  J. 

C 

JL  UO  1 = VRNKP1,  SC-' I 

V(I,J)  = V(-,J)  - PFOD*QPV  {?,:i 

120  CONTINUE 

C 

C END  LCCP  CVEF  T Hf  'c  LUMPS  OF  V. 

C 


no  r-.  non  non  non 


130 


CO  N T I N U E 


FNP  ICJP  CVEP  THP  “'“UPSHCLDEP  TP  A NS  PCP  HAT  IC  ^S  . 
200  CONTTNJE 


PFPrOTE  THE  SOWS  F V,  ACC  P'^ING  TO  THF  INTEFCHANGFS  CARFIED  OUT  OV 
“^PE  CCIII.'INS  CF  CSV  DUPIN';  XT'*  PECUC^ICN  PPCF  THE  EEFF. 


2 10  2 30  = 1,  NCCL 

CALL  UNSCFrf  1,  NCjL,  IPER**,  V ( 1 , J)  , TEFP  ) 
220  1 = 1 , NCOL 

V (I,  J)  = T £sr  C) 

2 20  CTNTINJE 

2 30  C'^NTTNUE 
RFTOPN 
END 
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6.2.  Svibroutiae  FRMORT 


i 

I 

I 

i* 


f 

j 

i 


I 


6.2.1.  Purpose 

The  subroutine  FRMORT  forms  the  explicit  orthoRonal  matrix  Q 

that  is  the  product  of  r Householder  transformations  applied  on  the 

left  by  the  subroutine  HREDL  to  reduce  an  m by  n matrix  of  rank 

T 

r to  upper  trapezoidal  form;  al ternativelv,  the  matrix  Q may  be 
formed . 

T 

The  first  r rows  of  Q (columns  of  Q ) form  an  orthogonal 

basis  for  the  subspacc  of  vectors  spanned  by  the  columns  of  the  original 

T 

matrix,  and  the  last  (n  - r)  rows  of  Q (columns  of  Q ) form  an 
orthogonal  basis  for  the  set  of  vectors  orthogonal  to  these  columns. 

6.2.2.  nescription  of  method 

The  matrix  0 is  the  product  of  the  r tnnsfonr.3tions  a.c 
applied  during  the  reduction  from  the  loft; 


Q = tl  H , 
r r-1 


H,Hi  , 


where  the  vector  cor rc.-spond ing  to  is  zero  in  components  1 through 

T 

j - 1.  The  matrix  Q is  the  product  of  these  transformations  in 
reverse  order: 


Q = H,H„  • • . H 

12  r-1  r 


A5 


I 


To  form  Q or  the  transformations  are  multiplied  together, 

starting  with  H in  explicit  form  — for  example,  in  forming  Q 
each  successive  transformation  is  applied  on  the  right,  leading  to 
the  recursive  definition: 


Q 


r 


H 

r 


for  k = r - 1 step  -1  until  1; 


Tlien : 


Q ^ q. 


In  applying 
Cure,  and  the  known 
(e.g.,  contains 

upper  left  corner)  . 
the  code. 


H|^,  full  advantage  is  taken  of  its  special  struc- 
form  of  the  partial  product  to  which  it  is  applied 
a (k  - 1)  by  (k  - 1)  Identity  matrix  in  its 
Complete  details  are  given  through  comments  in 


6.2.3.  Keywords 

Complete  orthogonal  factorization;  QR  factorization. 


6.2.4.  Source  language 

Fortran.  The  code  in  FRMCRT  has  been  checked  by  the  PFORT 
verifier,  and  is  WATFIV-compatible . All  variables  and  functions  are 
explicitly  declared. 
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I 

I 


I 

I 

[ 

i 


I 
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6.2.5.  Specification  and  parameters 
See  accompanying  listing. 

6.2.6.  Krror  indicators 

See  accompanying  listing  (the  description  of  the  parameter 
LERROR) . 

6.2.7.  An.xlllary  routines 

KR.MORT  calls  the  standard  functions  DABS  and  l.\BS . 

6.2.8  Program  si/e 

89  Fortran  source  statements. 

6.2.9.  Array  storage 

No  locally  declared  arrays. 

6.2.10.  Timing 

The  number  of  arithmetic  operations  required  to  form  Q is 

2 3 

of  approximate  order  2mr(m  - r)  + -^  r . 

6.2.11.  Further  Comments 
None . 


5 


47 


n n o 


S’lLPCIiriNt  FPrORT  ( NfiCW,  NPNK,  NDII,  QFV,  HVECl,  KUIOED, 

1 SCDIM,  CPTH,  lEFFCP) 

r* 

tNTFr.EP  NF-W,  N[^^K,  NDI",  f*OLOHr,  NODT**.,  IFPFCB 
DO'!"LE  PhtCISiCN  0PV(‘;ni‘*,  NPNK)  , HVECL  (NRNK),  0PTH(N0D:*!,  NHOW) 
C 
t 

f 

c 

C THE  JIJ380UTiNE  FRr'OFT  IS  USFD  (IN  CON.TONCTICN  WITH  THS 

C SHHPOHTTNr;  HBEDL)  TO  Fl*:,''  THE  EXPLICIT  CSTHCCCNAL  EATPIX  THAT  IS 
C TH’’  PFTPUCT  OF  A SPECIAL  Sr'OUENCF  OF  HOUSEHCIDFF  TP  ANS  FOP,  H ATTONS , 
r ''HLTiri.TE 0 IN  hITHFF  F'FWAPD  OR  PACKNARC  OFftP.  THE  S^OHENCE  OP 
C TPA  NSFOPF.ATI  J.NS  ARISES  WHFN  THE  SHBPGHTINF  HRECL  IS  USFD  TO  REDUCE 
C THE  'IMPIX  ghV  D UFFFP  TPAFEZOTDAL  FOP".  THIS  FPDUCTTON  INVOLVES 
C A SEO'iEVCC  OF  HCJO^MCiDEP  " I A N «:fC  P*1  A TI 0 NS  AFFIXED  CN  THE  LEFT  TO 

C URV  (FOSSIbLY  WITH  ■!'T  S CUUPNS  T NTEPCM  A «• ’.EF)  , ANT  «aY  WRITTEN 

C FtNRNK)  • * ’■(1)  ♦ (ntRr.UlEP  QPV^  = (UPPER  TRAPEZOID). 

C 

C EACH  F (0)  IS  A Hi.US£HCLrER  ''ATFIX,  OF  THE  FCRP 

C 

C I - ( MO)  * U(J)  TFAE  SP'S'-)  / BBTA(J), 

r 

C WHERE  y (J)  IS  A VECTCF  WHORE  FIPST  (J-1)  COP^CNENTS  APE  ZERO,  AND 
C 'STA(J)  IS  H SCALAS  NCR  "A  LI  7 1 NO,  FACTCP.  THE  SUPFOUTINE  HRSDL 
C CAPRIFS  OJI  IHIS  HFDUCxDN,  AND  STORFS  INFC  FH  ATTCN  AbOtlT  THE 
C TBA  NSF''F  .Y.M  ICNS  IN  CO''FACT  FOP'’  IF  THE  NATPIX  QPV  AND  IN  SEVERAL 

C AUyiT'^A'rY  VECTi,RS  (SFK  ""HP  D CCU  HE  NT  ATIO  N OF  HREDL  FO*?  DETAILS). 

C 

C THE  SUHRCUTIN"  pFF.''RT  USES  THE  OUTPUT  CP  HRECL  TC  GENERATE 

C ITTHEP  THE  MATRIX  C OZFI’^EL  AS  THE  PRODUCT  0”  THE  T PA  NSFO  RE  ATI  ON  S 

C IN  Fo:  WARD  .SD'B  , I.  S.  , 

/** 

C C = p (NPNK)  * . . . * P (2)  * P ( 1 ) , 

C 

C OP  THE  TK.ANSPObE  Or  0 (OT ) , WHICH  IS  GIVEN  ?Y  THE  PPODUCT  OF  THE 
C TEA  NSFCFMATIO.S  j IN  R V FP  S E C F DEF  , I.F., 

C 

C OT  = F ( 1)  • F (2)  • . . . • P (NRN  K)  . 

C 

C 

c THE  ao«s  CP  g (''OLU''Nr  of  o’:)  provide  useful  infcp>'ation 

C ABOUT  TME  ^uLUFN  ..FACF  F THE  ORIGINAL  , “ATP  IX.  lET  'NRVK'  PF  THE 
C r.STI“''T''D  NUriERIIAL  F .A  N K CF  TH'E  TPIGINAL  MAT’’IX,  AS  DE'’’EREINED 
during  EXECUriUN  t E HREDL.  THE  ‘•'TRST  NFNK  FJWS  OF  0 (COLUMNS  OP 
OT)  P'PM  rtN  C?THOGri,/>i  FASIS  F"' P THE  SPACE  SPANNEF  BY  THE  COLUPM.S 
OF  THE  OFIGINAL  lAIRIX,  AND  "HE  REMAININO,  ROWS  0 (COLUMNS  CF 
OT)  FOP"  AN  .RTH^GiNAl  BASIS  F ' F.  THS  SPACE  OF  V’^CTORS  ORTHOGONAL 
C THP  ''JLUMNS  CF  ^ riGINAL  MATRIX. 
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n n r>  n f'-  n n n n f . n r\  n n 1 r>  n i o <"1  rj  n o < < r.  o o n o o n o o n n o o o n r,  n o o o 


I 

1 


♦*  CAI'.UN  «•  Arvp--Ai.;i  IS  TAKEN  OF  THE  SPECIAL  FOPE  j?  THE 
TPANSFOPr.A'I  iCNf.  IN  ;-  Ff*l'*r-  "-HE  DpooOC",  SO  THAT  firOFT  IS  S’DT 
SHI  TABLE  FCF  FU  I " 3 I L Y J ’.’0  A OENFPAL  SFQHENCE  O'  HC'ISEHOLDER 
TPA  NSFC  Af.ATIOrJS. 


N FO  W - - 

INj-’EGEP,  INPHT  ONLY. 

THE  NUrHES  OF  CF  THE  OPIOINAI  '•ATRIX  C FV  THAT  WAS 

FFDICtO  BY  IF’i  SO’t)  UC  HI  IN  E MRFOL.  '«JST  3F  .01.  0. 

N R N K - - 

TSTEGES,  INIUr  C.  N L Y . 

THc  NC^BEP  . FAE-F.  ATIOKS  APPIIET  BY  HREDL  (7HF, 

ESri.NATFD  lANK).  "'HF  VALUE  OF  NPNK  FUST  NO"^  PE  ALTEREP 
AFlEb  £./  10  HP  Fl'l  AND  PSFCFE  EN'^FY  TC  FF-ORT. 

N c:  f1  — 

I'T'SCtP.  , INP  JI  CN'-Y, 

THE  DiioLAfED  S N PIF.FNMON  OF  THE  XATFIX  CFV  IN  THE  CAl.LISO 
S'JIi- FFCO  lA'' . .'FfT  PE  . G”.  NPOW. 

OFV  -- 

D'USLL  PREOISI""-  ARPAY,  0 OECLAPEO  P’‘HENSI0N  NEFF  PY  NPNK, 
AND  CoS  EEIUAL  DIHEK^'ICS  NRCW  BY  NPNK.  INIUT  <'NLY. 

THE  ''AfPIX  ORV  C pnESFCNDS  TO  THE  ORIGISAI  ^ATPIX  THAT  WAS 
F 'DUO  ED  BY  hPLDI.  TT*^  CONTENTS  '^CST  XDT  PE  ALTFPED  '.'‘"'FP 
E/IC  F50>'  HO  EDO  AND  '=tF''rE  F.NTtX  TO  '•RHOPn.  a tESCRIF'^IOS 
IF  THE  CO  NT-  NT.'  .F  QSV  IS  GIVEN  IN  H o DCCU ''EN'^  ATION  P 
POEDL. 

HVFCL  -- 

DOUoLE  rFECldl-N  VECTCE  OF  LENGTH  NRN»,  I MU'^  CNLT. 

T?:S  VEOT'  P MVECL  IE  GFNFPATFD  BY  fOEDl,  \ NO  ITS  CONTENTS  HUST 
MOV  HE  AMiPiD  A'’‘'SP  '.YIT  OPOH  HPEDL  A N F PFEOPr  ENTRY  TC 
FP-'orT.  A DtSCFir-TICN  OF  THF  CONTENTS  CF  HVECL  IS  GIVEN 
IN  fiii  DCru.X''NIA':ON  FCF  HPEDL. 

FULOFD  — 

IN^EOEP,  INPUT  CM  Y. 

TH';  SIGN  "ULORn  ISTICATES  THE  CROF.B  IN  WHICH  THE  TFAN-'- 
fORMAMCM  A .=  " T uF  H UI 1 1 X ■ D . IF  “UIOPO  =:  ♦ 1 , "HEY  APE 
.r’'L.  iPL  iFL  IN  r 'WARD  OFOEP,  I C IFOrUCP  Q.  IF  HULOPD  = -1, 
THEY  ARE  YULMPLTED  IN  R''VfR.SF  ORPEP,  TC  fRCD'JCE  Q IPANSPOSE. 

NODI.''  -- 

INTEGEo,  INTOT  ON'!  Y. 

THE  Dt-C'tARLD  «'''•«  CIH’N.iI5S  OP  THE  EA'!*FIX  CFTH  IN  THE  CALLING 
C S''H-FP  'GRA  .''U.'T  P:  .GS.  NPOW. 


i 

I 


o o o o o n n o n n fi 


C 

C CFTH 


r DOUBLE  PpBJISnN  .SF  = AY,  D3CLAP»')  CIM’'N’.SION  NODIB  3Y  NPOW, 

C A*!D  CONCEFl'JAL  T I.'!  EV "C  V V'^'C'W  BY  N'’OW  . OHIPOT  ONLY. 

C ■;%’  EXIT  FFC.Y  FFY.r?':,  THE  ''A'^PIX  ORTH  CON^'ATFS  FITHF?  Q 0”  OT 

C (0  rF.ASSFCSE»,  DtFFNPTN^  ^ ‘i  THE  S lOK  0’^  rriLCPD. 


C 

C 

C 

C 


i. 


r 


IFSFTP  -- 

INTEGEP,  CUTFIII  CNLY. 

FBPOF  INDIJAICF,  5”T  PDFIVC  EX’='CtITION  0’=’  FFYCPT,  WITH  THE 
FOLIPWING  ^OSSIBl'^'  VAItlFG. 

IE  PROP  ^ 0 : '"'=FC''’S,  T'=’ PHIN  ATION  . 

LEfRCP  = 1 : INVALir  T\-pnT  FA'’A’1?rFP. 

LF.FPOP  - 2 : vfXK  IS  7.  F?'^ . IF  LEFP0P  = 2,  0 IS  SET  TO 
■"HE  in2N'’*TY  HATPIX. 


A'JTHDPS; 


HiP  JAFF"'  H.  WF:g«T,  STEVEN  C.  GLASEY.AN 
SYSTIYS  OFTir  I7A 'TION  LABCPATCFY 
DTFAFTKEN’T  CF  CPEPATIONS  PFSEARCK 
SIANFDFD  UVIVE  P'^I'^Y 
SIAN  FTP  r,  CALI’='^PNTA  9tt3C5 


«**  p^stE: 


DECFF.'^EF  1977 


C 

PECLAPATI7N  OF  LOCAL  V A P"' AP  L'=‘ S . 

INTSGS3  1.  J,  JP1,  K,  KK,  KP1,  NPNEP1,  KPFKM,  NF^WHI 
CCL'SIE  PPECISICK  F'^TA,  FAC^P,  PCTIFD,  FCTJ,  FCTPNK 

DECLAEAT  ION  Cr  SIANDA"P  FfNCT'^ONS. 

IKTEGE3  IA5: 

D'^'JELE  r EEC  IS  ION  CABS 


TEST  FOP  EFEOR  IN'  IN’FI!!  ?APf''ETFFS. 

L.EE'^'ip  = 1 

(NFOW  ) .OF.  N’’NF  .LI".  ^ .OP.  N^^ir  .IT.  N'^rw  .OR. 

1 NODI*  .IE.  C ."'E.  NFNK  .GT.  NFOW  .OF. 

2 I A3S  ('IfJLOPOl  . NF.  11  '’"'T'lFN 


INITIALIZE  ORTH  TO  THE  '"'W  HY  N F'^'W  IDENTITY  NATFIV. 
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DO  20  J = 1,  .iPOW 

00  10  r = 1,  «Fow 
-'F'i.  1(1 . j)  - c.oc+o 

10  CCNriN'IE 

OPTli(J,J)  = 1.0D+') 

20  CCNTTNUE 
IPFBOF  = 2 

IF  (NENK  .SO.  0)  FFT'JPN 


GENEPATP  ? (NHNK)  E.tfLICirLY,  SIMCF  ITS  FOR.*;  IS  THAT  CF  A SINGLE 
HOUSEHOLDS^  ih  A t.3  F "lE  HA  T ' CN. 


IF  (HVEOL(NPNK)  .FO.  O.OP  + 0)  GO  TO  HO 

9ETA  = D AES  (H  V£CL  (NFMK)  1 

FACTOr  = 1.CD*C/?ETA 

OFTF  (N.-vN  E.  NR.VKi  = I.Oi+C  - B’='TA 

IF  (NPN'K  .20.  N?OW)  GO  '~C  00 


.'•ILL  IN  Trih  LOWE?  RIGK"'-HAKD  (NRCW-NFNE.)  PY  (;4riOH-NPNF)  SQUAPF  OF 
P (N'FNKI  ; THE  EE.HAINEiF  C F PfKFVK)  IS  SIFPIY  THE  ICFVTTTY. 


TMF  NRNK-IH  ROW  AND  COLUf-N  CF  P(NFNK)  APE  FCP.^KD  SEEAFATELY 
SINCE  THE  NRNK-TH  CO  HPO  VT  rp  fl(NENK)  I.S  .STOFFD  IV  PVECL(NRNK). 
THF  FFIAINING  COYP'^NFNTS  CF  U(NRNK1  A P STCp-^O  5FI.0W  THE  DIAGONAL 
I’.'  TH'’  NRSK-r.H  CZL'JYN  C F CRV. 


NFNKn  = NRKK  1 

FCTPNK  = FACT  IP  •HVECL  (NPNK) 

DC  3C  I = NFNKPl,  Vie  w 

OFTil  (I  , NFNK)  = -r  CTFNK*0'’ V (I,NFNK) 
CF^!^(NR^^,X)  = ■'FTH  (.',  NPNK) 

30  CONT'NIJE 

IF  (NRVKEl  .g:.  NFCWI  GC  Tj  60 
NRCW*'  1 = NRCW  - 1 


THE  ljjP  siAiErEr.;  50  conrurrs  all  elfisntj  of  the  sub-hatrit 
except  the  NFOW-IH  DI  AGON  A I,  ElEMCNT. 


DO  50  J = NFVKP  1,  NFCW"! 


n n r n <->00  o o <')  n n t d n ct  n n n n f i r 000 


JP1  = J + 1 

FCIJ  = PACfOF^'wFV  f.I,  NFNK) 

Cr>  4 0 I = JP1,  NF._H 

= -FC:J*0K  V {I  ,NFNK) 

UO  CON  TIN 'IE 

CHr;i(J,J)  = 1.0C  + '^  - FACT ’>P’‘ORV  (.1  ,NPNK)  **2 
50  CONTTN'JE 

DEFINE  THE  NFOi<-TH  Di.AGCKAL  ELFr.HNT 

60  OPFH  (N.<  CW,  NFtW)  = 1.0+0  - ' A CTOP*  QR  V ( N P CW  , A F A K ) * *2 
70  CONTINJF 
HO  LEFBCP  = 0 

IF  (NRNK  ,SC.  1)  FLTURN 


THF  l''OP  TO  SIATEME.VT  200  HdLIPLIES  TOGETHEF  THE  FEHAIHING 
H0USEH01.D£,P  rRANGEcF'-ATIONS.  BEFORE  BEGINNING  THF  LOOP,  ORTH 
cnN'''ATVS  P(NI-.NK).  EACH  '’I';F  THROUGH  THF  LCOP,  TH"'  CURRENT 
OPTHOGCSAL  rAIFIT  iJ  ''.ULTirTISD  EITHER  CN  THE  BIGHT  (IF  HOLOPD  = 
OR  ON  THE  LEFT  (IF  .'lULOrU  = - H 3Y  A HOUFEHCIDER  T FANSFOPM  ATION 
CCPRF5P0NDING  TC  A VECTOR  WI'TH  OAE  .-ORE  NJN-TERO  COMP'^VENT. 


♦ 1) 


NPNKTI  = NFNK  - 1 

DO  200  KK  = 1,  NFNFH1 

FOB  K = NFNK-I  STEP  -1  UNTIL  1 

K = NFNKn 1 - KF  + 1 

TEST  WHETti-:?  THF  K-TH  '^P  Ml  S FORM  AT  ION  WAS  SKIPPED. 

IF  {HVECL(K)  .iO.  '^.OD  + Ol  GO  TO  200 
PET  A = LAF  E f HVEOI  (K)  ) 

Fi'.CTOR  = l.CD+^/Rt:A 
IP  (M'lLEFr  .IT.  J)  GO  TO  150 


THE  K-TM  T;  A NEFCF- ATIC  N IS  TO  EE  APP1I'’P  CA  THE  RIGHT. 

C 

C ■'•HE  FIRST  (K-ll  FOW£  rf  pHF  CUFRpAT  C AFE  CNALTEPFD 

C WMEN  P(K)  S Ain:  ED  FFOM  ^ HE  FIGHT,  SINCE 

C C)(I)»*T»U(K)  = 0,  I = 1,...,K-1,  WHERE  0(T1**T  IS  THE 

O'  ■'■-Til  FOW  ' f'  TH::  C'JFFLf.T  0.  THE  FIRST  (K-1)  COLUMNS  OF  THF 

C E'lISTING  0 ARc  ALS  UNALTEFFD,  SINCE  THE  FIRST  (N-1) 
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C Cf"lt'3N£N‘r.-i  CF  !)(!')  APF  ZKPC,  AKU  'THE  FfPST  K ROWS  J? 

C THE  CliPPEiir  Q .J'.PT  PCWS  !>F  THE  I^'F^TiTY  .-AT  FIT. 

C 

C CCTIpnTe  THJ  K-TM  POW  ry  ‘UDIFIEP  O.  which  is 

C THE  h-IH  PJW  or  P<K),  EE'-'AHSi:  TH  P K-TH  RO  W 0 ^ 

C THE  UN'!DDI?IFr.  0 IS  THE  K-TH  5CW  THE  TDFNTriT. 

C 

C 

c 

OFXH(K,K)  = 1.CC  + 0 - FACTOF’HVhCE  (K>  ♦ *2 
KP  1 = K » 1 

PCTIND  - FACTC.'*HVECL  (K' 

TO  no  J = KP  1, 

JPrM!K»0)  = - FTTINr^O*- Y IJ,  K) 

no  COKTINiJZ 

L 

C THE  LOJP  VJ  STAIFnENT  nO  AL'''FF5  «='' W S TFPOUOH  Spf'W  'F' 

C THE  PPEVIOJSLY  COEFIITID  0. 

C 

Dj  no  : - Kfi,  NFcw 
C 

C FGhH  the  £C»i.4’^  TP' PUCT  ^F  THE  I-'T!  F '' W CF  THE  !IN  F.  JD I Ft  ED 

C V'  and  the  VrCTCP  OfK).  SINCE  THF  K-TH  CrTO-N  TP  TM ' 

C ,J  N w iOlFTE  0 g IS  ZErO  IN  FOWS  K ♦ 1 TH’OUOH  NROW,  THE  '^ERr 

r )?  THE  INNEP  PPOI'iJCT  C I PPESPTN  CT  N C THIS  CC'^ONENT  CF 

C OIK)  IE  nTITT'^P. 

C 

DcrppD  = o.''r*c 

00  no  J = KP1,  NPOW 

UCTHPr  = CCTPRO  ♦ ' PTH  (1,0)  *CPV  fJ,K) 

120  CCMINHE 

C 

C iODIFY  TH:  r-"H  r-CW  HY  SiJDTPACTTNG  an  appropriate  '■.I'LTTPLF 

C Or  IHE  H.:  CSr  HCLITP  VFCTOP,  UfK). 

». 

Fj;j  - rA  JT-c  ■‘OC'PPO 

ORrL(^,K)  - .■:»THf:,K)  - FCTJ  *K  VECL  (K) 

00  1 JO  J = KP 1 , NF'  W 

aRT'Ml,J)  ^ CPIH(T,.1)  - FCT  J*OF  V (J,  K) 

110  Jv  N T IN ')  1 

C 

C STATEMENT  14C  TJ.[)3  THF  LOOP  CVEP  THF  MTEFFC  ROWS  OF  g. 

C 

14  0 N r : T.  U E 

GO  ?C  200 
C 

ISO  contik'jf: 

r* 

C 

c 


S3 


o n ~ o n o rt  n r.  n i~,  r.  r.  o n o ri  o n rj  o o o n n c c r>  o n r;  o 


TI’E  K-ni  IfiA  NSrCR*’\T'CN  IS  TO  PE  APFLIED  ON  THE  LEFT. 

T“fc  FIPST  (K-1)  COLU/N^;  OF  THE  CHRRFNT  CT  ARE  ONALTEPED 
HPtN  P(K)  IS  AfTL’^Sr  FPuV  TH5  LEFT,  SINCE 

0T(*)  ’■*T*U  (K)  = : r 1 K-1,  WHFPE  OT  ( I)  IS  THE 

I-TH  COLH.-N  ''F  THt  C'FPFNT  OT.  THE  FIRST  (K-1)  ROWS  OF  THE 
EMSTING  Cl  AEE  ALSO  ’JNAI.TEPEO,  SINCE  THE  FIRST  (K-1) 
CCnPONENTS  uF  IMF)  ARE  7.KPC,  AND  IMF.  FIRST  K C0LII''NS  OF 
THE  C'JRRENI  01  A:  , v.."*!  H « S S OF  THE  IDENTITY  FATFIX. 

COIIPIITE  THE  K-IH  J'LUFV  07  THF  HOCIFIED  OT,  WHICH  IS  SIHPLY 
THE  K-IH  CjLUEN  OF  P (F  ^ , BKCAIISE  THF  K-TI-  CCLUHN  OF 
THE  IJ  N.IODI  El  FD  01  "S  1 HE  K-TH  CCin.MN  CE  THF  treNTITY. 


CF;'H(K,K)  = I.CC*'^  - FA''T'>R*HVECL  (K)  **2  ] 

KP1  = K ♦ 1 

FCTIND  = F.ICrC  B«l:V  ECL  ('■') 
n 160  J = KP1,  NP"x 

0FTH(J,K)  ^ -i-CriNT'O  V(J,K) 

160  CONTIN'JE 

THE  LOOP  STATT.YENT  1^0  AITEFS  COIMHNS  K*1  THPOUGH  NRCH  OF 
THE  IREVICJSLY  CCEPl’T'F  O'". 

C?  1M0  J = KPI,  NROW 

EORE  TH  : SgAL.SR  PrCC'f’t.  2F  TH.E  CCIMEN  OF  THE  HNroriFIFO 

2T  AND  I'HE  VECT'r  ''(K)  . SINCE  IHF  K-TH  PCH  CF  THE 
UN'.IODIF  lED  OT  lo  2FF'''  IN  CCLM-NS  K*1  'T'HFOMGH  NPOW,  THE  TE  Pfi 
)F  IHE  INNER  ri'C"  " "'H  ? ES  P 0 NDI  HG  TO  THIS  COMPONENT  OF 
U (K)  IS  OMITTEP. 

OOTPRD  = 0.00+'’ 

DO  170  I = KPI,  N'  . W 

DCTIRD  = DOIPRP  ♦ "P  T H ( I , J)  *0R  V ( I , K) 

170  CCNTINHE 

.MODIFY  IHE  J-IH  COItJ'''-  FiY  SUBTRACTING  AN  APPROPRIATE 
MULIIfl.-;  CF  TH-'  H''USE.HOLDER  VFCTOP,  U(K). 

PCT.I  = .‘•’ACT.IR  *DOTrFD 

ut.H(K,J)  = rTH(i<,,1)  - 7CTJ*M  VECL  (K) 

DC  1 dO  ' = K F 1 , N r ' W 

7PTH(I,J)  = '.R"H  (!,.])  - F CTO*OP  V (I,K) 

IHO  -'(.NIINUZ 

STAIErtENT  ImC  ENDS  iHI  L-'OP  .OVEP  THE  A ITEREO  COLUMNS  ""F  OT. 


6.3.  Subroutine  HMULTC 


6.3.1.  Purpose 

Tlte  subroutine  HMULTC  (for  Householder  transformation  multiplied 
over  £olumns)  constructs  a single  Householder  transform.atlon  to  annihi- 
late a block  of  components  in  a particular  column  of  a given  matrix, 
and  applies  tlie  transformation  to  the  remaining  colvimns.  The  process 
of  defining  the  transformation  is  described  in  Section  2.1. 

6.3.2.  Description  of  method 
See  Section  2.1. 

6.3.3.  Keywords 

Orthogonal  transformation;  Householder  transformation. 

6. 3.  A.  Source  language 

Fortran.  Ttie  code  in  HMULTC  has  been  checked  by  the  PFORT 
verifier,  and  is  WATFIV-compatlble . All  variables  and  functions  are 
explicitly  declared. 

6.3.5.  Specification  and  parameters 
See  accompanying  listing. 

6.3.6.  Error  Indicators 

See  accompanying  listing  (the  description  of  the  parameter 
LERROR) . 
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J 


r 


! 


j 


I 


\ 

I 


% 


1 

1 


6.3.7.  Auxiliary  routines 

HMULTC  requires  die  standard  functions  DABS  and  DSQRT. 

6.3.8.  Program  size 

37  Fortran  source  statements. 

6.3.9.  Array  storage 

No  locally  declared  arrays. 

6.3.10.  Timing 

riie  number  of  arithmi'tic  operations  necessary  to  con.struct 
a single  normalizeii  Householder  transformation  to  annihilate  k 
components  of  an  r>-vector,  and  then  transform  2,  other  m-vectors 
includes  k + ?.  divisions,  and  of  order  2Z(k  + 1)  multiplications/ 
adili  tlons . 


6.3.11.  Further  comments 

Representation  of  Householder  transformation . A Householder 
matrix  is  defined  in  terms  of  the  non-zero  vector  u,  by: 


H(u) 


where  p (tfie  scaling  factor  for  tlie  transformation)  = 1/2  Hull*’.  To 
represent  such  a transformation,  it  might  appear  that  only  the  vector 
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u needs  to  be  retained;  Itowever,  it  would  usually  be  considered 

2 

wasteful  to  recompute  1/2  SuN  several  times,  and  thus  one  might 
wish  to  store  the  quantity  B as  well. 


The  transformations  used  in  reducing  a m.atrix  to  upper  trapezoidal 
form  from  the  left  liave  a special  structure  (see  Section  2.1);  in 
particular,  each  transformation  may  be  regarded  as  reducing  the  "first" 
column  of  a remaining  matrix.  Consequently,  we  need  only  consider 
representing  the  Householder  transformation  that  reduces  a single 
vector  y with  Z components  to  a multiple  of  e^^,  l.e.,  that  annihi- 
lates all  but  the  first  component  of  y: 

D 
0 
0 

0 

Because  Euclidean  length  is  preserved  by  a Householder  transfer-  ^ 

2 2 I 

mation,  p = Nyll.^.  Multiplying  out  the  first  expression  in  (17),  it 

can  be  seen  that  the  Householder  vector  u is  a multiple  of  (y-pe^^),  I 


i .e. , 


- P] 


for  some  scalar  y.  To  avoid  cancellation  In  computing  n,  the  sign 
of  p is  always  chosen  to  be  opposite  to  that  of  yj,  so  that 
p = -sign  (yj^)llyll,  where  slgn(y^)  = 1 if  y^^  2.  other- 


wise. 


Tlie  Householder  vector  u can  be  stored  very  efficiently. 

First,  since  components  2 through  i of  v will  be  annhilated  by 
the  transformation,  components  2 through  f.  of  the  Householder 
vector  can  be  stored  in  these  components  of  y;  hence,  only  one  addi- 
tional location  is  required  to  store  the  first  component  of  u. 

Second,  the  scalar  y In  (18)  can  be  chosen  so  that  the  scaling  factor 
d need  not  be  stored  separately  (Stewart,  1977).  Tf  y “ tne 


value  of  3 is: 


S = i llull^  = 


^ (y^  + 2iy^i  lly'l  + Hyll^  + y?  + •'*  >’}) 


llyll  + ly^l  Nyll 


= 1 + 


If  Y - then 


= Yj/llyli  + sign  (y^),  so  that 

loil  = 1 + ly^l/llyll  = 6. 

Thus,  the  Householder  vector  u is  given  by: 


u,  = + sign  (y, ) 

i N II  t 


y. 


(stored  separately) 


u = — ^ , j = 2,  ....  2.  (stored  in  j-th  component  of  y) 

j llyli 


This  scaling  of  the  Householder  vector  eliminates  the  need  to 
store  8 separately  or  recompute  6.  Furthermore,  the  Householder 
vector  is  normalized  so  that  its  length  lies  between  and  2; 

this  latter  property  is  useful  in  avoiding  overflow  or  underflow  when 
applying  tlie  transformation  to  other  vectors. 
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StHFL'UrisF.  H-nULlC  ( IPCrf,  ICOL,  MPOW,  MCOL,  NDTK,  QRV, 
1 CNiri,  U.MPti',  LERTCP  ) 

C 

IN-'i'SER  IF>.  W,  NF';^,  NCOL,  NDTM,  L-RPOF 

Drf'RLE  FPEi-ISIOS  CMOPM  , 'JMPHfi 

DuUBLt  PRECISION  (^PV(NDIM,  NCCL) 

C 

C 


TH=:  SUBRCUTINE  ‘^KULIC  (FOF  fcCUSEHOIDEF  T E A F3  FC  PE  A ITON 

r. ULTIPL :ei)  over  cciu-.nsi  wii.r  coNsrpucr  th’^  housfhcldef 

rS.'.NErOP'lA  ION  DEJIGNEP  ZC  ANNIHILATE  r''LEEFN''.S  TPOrf*-1  TriSOUuH 
NPi  W "F  CJL'IIN  IC^yL  r T !i"  rAT^^X  O'.-V,  A sn  TK3N  fPPI.V  THIS 
TPA  N5F0P  "A  I ! CN  TJ  'CLHEVo  ICCI.+  1 THPC’^fi  NCCl  CF  CPV. 

TK-  COUPON"  NT  •.  F ''ME  VECTjP  PF  FIN  TNG  THF 

'FlHSEF'C!  Lt,r.  :i  AN  ;?.PNAT  IjN  15  ST1FFJ  IN  (lALFHA.  I"''  VECTOR 
I?  NOPNAL'.  bi  TIU:  THE  F.AGNIT'JDE  OF  L'AIPHA  IF  MS.''  '^“K 
rCAt.’’NC  FhCT  j&  I :l  THF.  TFANSFCPX  ATION  (EEE  the  SX'^ESNAI. 
D''Cr)F:F.NTA'ir  N FO?  FUFTHEP  OE'T'ATLS).  C.ONPONFNrs  IP'^W+1 
rhPCrOH  Nrc-  >,F  THE  VECTOF  APE  STORED  IN  THE  C C F F Ef  TO  NPT  N G 
tLEXENTE  Jr  OCLlI-'t.  IC'"!  .'<F  C?V. 


I 


( 
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THE  F'PfAL  FABAE-TEE:  .PF  FlXfir^c  AFE: 


INFEGEF,  INFHI  'NLY, 

1FC«  IS  THE  TNOF-X  FOP  '''Ft  STAFTINC  EIFrEM 
liC'JOEFiCLDEF  VECTCP.  EHST  PE  .GF.  1. 


INTEGSf,  INPilT  ''KLY. 

:C01  IS  y ME  INPEX  CF  THE  CJLHMN  MfON  WIMCH 
:F  ANSFORMA' 'CN  *S  RASED  . XUS''  3E  .CF.  1. 


CF  "HE 


THE  HOHSFHOI.DFP 


NF  .W 

t IN'!  EGER , INtt'r  .'NLY. 

(•  :<FO«  IS  th::  infev  (F  the  last  element  ir;  the  cdlhhns  of  osv 

C rc  BE  AFFFCTEl,'  - Y THE  H''MS'"H  CITE  F T F AN  EFCRX  ATION  . 

r XU  SI  BE  ..AK.  It'.  Vi. 

C 

C N CC.  L 

C iNTHGEF,  IRptJ-  K I Y . 

C SCCL  'S  THE  INPEX  t F THF  LAST  COIHFS  TC  ViHICH  THE  HG  OSF- 

C HCIDEP  . PANS'-.  ^XATItV  IS  APPLI'O.  NHST  HF  .GE.  ICOL. 

' C 

C NOTE  

--  INTEGER,  IKPUl  OMY. 

C NDIX  IS  I Mb  '-XIFFNALLY  CFCLAR»^C  POW  DIFENSTCN  np  THF 


I 
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on  nnnoonoooooooorinonrionoonooir'ooooooorioooononoon 


r 


I 

I 

<a:fix  of.v.  .t;r.  nfcw. 


ORV  — 

UDUBLE  PR  Pl'I-SICA  ARRAY,  OK  CONCEPTUAL  DIMENSION  NPOW 
i3Y  NCOL,  AN''  PtCLAFPD  DIMENSION  NPIM  EY  NCOL. 

INEUT  AND  OUTFOT. 

ORV  IS  iHP  MArFIX  10  HR  TRANSFORMED.  CN  EXIT,  ROWS 

IPuM  THROUOll  NROW  f F COLUMNS  I CO  I THROUGH  NCOL  OF  QF.V  WILL 

HAVE  bEEN  ALTERED. 

CNCRM  

DOUBLE  PR  EC  IS  ION,  INPUT  ONLY. 

CNUHM  IS  IhE  '^"01  IDEAN  LENGTH  OF  THE  FErOCEC  COLUMN  ICOL 
.<:-’C.I  NM  NO  WITH  ’^'LEMRN'"  IROV)  A.  N C ENDIKO,  WITH  ELEMENT  NPOW. 
■MUar  BE  .GT.  o.or+p. 


UA^  FH  A 

DCUJLE  EFECCICN,  CUTPHT  ONLY. 

UAIFHA  JCNr\:<S  .HF  IP  W-TH  ELEMENT  OF  THE  FiCUSEHOLDEP 
VFCTJF.  TH'=:  ABSCIUTE  VALUE  OF  UALPHA  IS  ALSO  ’’HE 
SCALING  FACT'  F FF  THE  TR  A NSF PMA  T I ’N  . IF  UALPHA  IS  EXACTLY 
ALRC,  TIU  TF  ANEFORM  AIION  WAS  SEIPFEB. 

LEPPOR  

iMEGLF.  MlTFUl  (NI.  Y. 

L"'!  FCP  IS  run  E?R  F FLAG  AND  HAS  TWC  E^SSIBLE  VALUES. 

LERPOF  = C ; ’IOF''AL  FETUFN,  NO  ERRORS  FOUND. 

LEFbJR  = 1 ; ='pr-'fJ  IK  INPUT  FAFAMETEFS. 


AUTHORo:  LAFGAR^T  H.  WFIGHt,  STEVEN  C.  GLAE'MAN 

SYSTLME  P'il  M:TA''’I0N  LABCRATCRY 
DEbAF’'M.EK  " Of  OPERATIONS  RESEARCH 
S”ANFiRD  UNIVrRST'''Y 
STASF'PD,  CAII''^RN:A  9U3C5 

***  CATE:  DECEMPEF  1977 


DECLARATION  CF  LOCAL  VAFIAPLFS. 

INTEGER  I,  ICLFI,  IRWM,  KCOL,  KROW 

double  PRECIS  ION  SETA,  DOTPRD,  ELEM,  FSIGN,  GAMMA 

DECLAfATICN  C'  STANLAFD  FUNCTIONS. 
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AFFLY  I’Hc:  rFA.'«5F  S.iATlCK  T.  TH!=:  PEYATNING  CCLCr.NS  OF  QR'l.  IF 
iCCl  li  IHF  CNLY  'JCLUHM  '0  " E TF  A NS  FCP  **.  E P , EXIT. 


I 
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IF  (ICOL  . EQ.  NCCL)  FETUFN 
I Cl PI  = ICOL  ♦ 1 

BETA  = DAhS(UALPHA) 

<4  30  KCCL  = IJLF1,  NCOL 
C 

L,  F'^ND  TFE  DJ;  PRODUCT  fDOTf’PD)  CF  IH  E HOUSFHOLDEP  VECTOR  AND 

C ■r’HE  COLUMN  1.)  RE  ! . A N S F'' R r ED , KCCl.  ELE.IFNl  TPOh  OF  THE 

C H.'USLHCLDLR  VL'l.F  IS  ST'^’En  IN  UALPHA  AFC  SC  IS  TREATED 

C 5*‘PAHArEIY. 

C 

DCrpFD  = UALPHA  • OR  V f TP  ) W , KCC  L) 

DO  410  KFc«  = " FWPl , NROW 

DOT FRD  ^ DClPPn  + OR  V (KRC W , ICC  I ) • QF V (K PCW, KCOL) 

410  CCNIiSUE 
C 



r 

C APPLY  THE  H0U3EH0LDEF  P.A  NSFOPH  ATION  TO  THE  COLUMN  KCOL.  THIS 

C IS  EOUIVAIEN'T  To  S UBTF  ACTIN’C,  OUT  A .MULTIPLE  OF  THE  HOUSEHOLDER 

C VECTOR  FFCM  IHE  COLUMN.  jNLY  ElEFENTS  LFCW  THFCUGH  NRCW  OF 

C THE  CClUFN  APE  F.<PLICITLY  MCDIFIEC,  SINCE  THE  FIRST  (IFCH-1) 

C C'^.MPCNENIS  CF  THE  P'USEH  ' LDEH  VECTOR  A FE  ZERO.  THE  IROW-TH 

C"HPCNEST  IS  TR’^ALTC  SEPARATELY,  SINCE  THE  IROH-TH  COMPONENT 
C O'"  i'rf..  IiLU:EHCLD£F  VEC10=  IS  STORED  TV  UAIPHA. 

C 

C 

c 

GAMMA  = DC.PRD/bFTA 

OP V( IROW  ,KCO  1)  = OF  V a RON , KCCL)  - GAMMA  • UALPHA 
DC  4:^0  KRO^  = .RfcFI,  NPjW 

OR  V (KF''W,  KC'L)  = OF  V (KFOR  ,KCOL)  - GAMMA  * C FV  (K  RO  M , 1 C'' L ) 

42''  CCNIINUE 

STA''’£.'1ENi  4)0  cNTS  THE  LOOP  ''VER  COLUMNS  ICOL«-1  THROUGH  NCOL. 

C 

430  CONTINUE 
FE'^UPN 
E ND 
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6. A.  Subroutine  UREDL 
6.4.1.  Purpose 

The  subroutine  HREDL  (lor  Householder  reduction  from  the  J[eft) 
is  designed  to  reduce  a general  real  matrix  to  upper  trapezoidal  form 
by  application  on  the  left  of  a special  sequence  of  Householder 
transformations  (the  process  is  described  in  detail  in  Section  2.1). 
The  principal  use  of  HREDL  will  be  as  the  first  step  of  solving  .a 
linear  least-squares  problem. 

6. A. 2.  Description  of  method 
See  Section  2.1. 

6. A. 3.  Kcjvords 

Householder  tran.s  format  ion;  Householder  reduction  to  upper 
trapezoidal  form. 

6. A. A.  Source  language 

Fortran.  The  code  in  HREDL  has  been  checked  by  the  P’^ORT 
v'erifier,  and  is  WATFIV-compa tible . The  current  code  has  been  written 
for  double  precision  floating-point  on  an  IBM  360/370;  the  machine- 
dependent  constant  EPSMCH  is  set  in  a DATA  statement,  and  should  be 
altered  to  correspond  with  tiie  given  computer.  Ail  variables  and 
functions  arc  explicitly  declared. 


I 

k 


65 


6.4.5.  Specification  and  parameters 
See  accompanying  listing. 

6.4.6.  Error  indicators 

See  accompanying  listing  (the  description  of  the  parameter 
LERROR) . Several  precautions  against  underflow  have  already  been 
included.  However,  in  extreme  cases  it  may  be  necessary  to  add  further 
safeguards;  comments  in  the  code  indicate  where  such  changes  might 
be  made . 

6.4.7.  Auxiliary  routines 

HREDL  calls  the  subroutines  H>R'LTC  and  LENSQ.  It  also  requires 
the  standard  functions  DSQRT  and  lABS. 

6.4.8.  Program  size 

137  Fortran  source  statements. 

6.4.9.  Array  storage 

No  locally  declared  arrays. 

6.4.10.  Timing 

The  number  of  arithmetic  operations  (additions/multiplications) 
necessary  to  reduce  a full-rank  m by  n matrix  (m  ^ n)  to  trape- 
zoidal form  using  Householder  transformations  is  of  approximate  order 


r 
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6 . A . ] 1 . Fur  t lie r comnien  cs 

Column  interchan^iij! . If  ■»  ■•■'U.'mn  interchaniie  strategy  is  used 
In  HREDL,  the  relei’ant  roliimns  arc  physically  Interchanged,  and  the 
exchange  is  recorded  iti  the  integer  array  IPERM,  whose  k-th  component 
givi'S  the  index  in  the  original  matrix  of  the  k— ch  reduced  column. 

For  example,  let  n 4,  and  Let  the  origin.-’l  matrix  A be  partitioned 
into  its  columns; 


A = [a, 


If  the  'PERM  array  is  (A  2 13).  the  matrix  actually  reduced  is: 


AP  - la. 


a,  a ] 

j.  j 


Definition  of  "size"  oi  a reingining  e.olutnn.  When  using  a 
column  Intel  change  strategy  -in  Cran.s  forming  from  the  left,  Mie 
"lai'gest"  remaining  columr.  i.s  chosen  to  be  reduced  at  each  stage. 
However,  it  is  impossi'fie  to  ilevise  a universally  applicable  defini“ 
tion  of  the  "size"  of  a remaining  column  (.see  Section  A. 2).  There- 
fore, three  column-based  options  are  provided  with  HREDL  to  allow 
some  flexibility  in  tiiis  deiiiiiMon.  The  size  of  a remaining  column 
can  be  specified  as: 


(a)  the  ratio  of  its  curient  length  to  the  original  length  of  the 
full  column.  This  measure  is  appropriate  when  the  elements  in 
each  column  have  a fairly  consistent  scaling  with  one  anotlier. 
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but  not  necessnrily  with  nil  cChor  columns  — for  example,  if 


each  column  represent-;  the  observed  values  of  a particular 
quantity,  but  the  scaling  of  each  quantity  may  vary. 

(b)  Its  current  length.  Tliis  is  the  most  obvious  definition  of  "size," 
and  is  appropriate  when  .1  uniform  scaling  exists  among  all 
columns  — for  example,  the  matrix  may  have  been  scaled  before 
calling  HRKDI.. 

(c)  The  ratio  of  its  current  lengtii  to  a user-provided  scaling  factor. 
This  choice  allows  some  correction  to  be  made  for  badly  scaled 
data.  It  should  he  emphasized  that  this  option  does  not  mean 
that  a weighted  least-squares  problem  is  solved,  but  only  that 

the  user  can  exercise  some  control  over  which  columns  are  considered 
"largest."  For  example,  this  option  could  he  used  to  encourage 
(or  discouragel  the  selection  of  p.irticular  columns  in  finding 
a linearly  independent  set.  or  in  the  early  stages  of  the  reduction. 

Definition  of  "ne^H^tlh’.;" . The  test  for  whether  a particular 
remaining  column  is  "negl  i I’.ih)  e"  is  of  the  form; 

(size  oi  column)  toler.ince  , 

and  thus  involves  the  spec  i t i*  at  ion  oi  a toler.mce.  Infortunately , tio 
universal  m.agic  value  of  such  .1  ti.li>ranco  exi.sts;  the  design  of  HRF.Dl. 
accordingly  includes  .1  procedure  for  defining  "negligible  that  is 
based  on  a combination  of  nser-provi iled  and  machine-dependent  informa- 
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Kirst,  the  parameter  COI.FPS  allows  t lie  user  to  control  the 

definition  of  "n«‘,’,l  iy.lhle, " based  on  known  properties  of  the  data 

and  on  the  nature  of  the  problem  to  be  solved. 

In  particular,  COLKl’S  should  reflect  the  accuracy  of  the  data 

fwith  respect  to  ihc  ciiosen  definition  of  size).  For  example,  if  the 

mairlx  consists  of  observed  values  with  at  most  A figures  of  relative 

-A 

accuracy,  COI.EPS  should  be  no  less  than  10 

The  type  of  problem  to  be  solveo  also  affects  the  choice  ot 
COI.EPS.  Rouj'.hly  speaking,  if  the  residual  vector  of  a least-squares 
problem  is  not  "small,''  the  conditioning  of  the  prtihlem  is  proportional 
to  the  squared  reciprocal  of  the  length  of  the  smallest  remaining 
column  (see  Stewiirt  . 1973,  for  a precise  statement  .nid  details). 

Thus,  for  numerical  stability  in  a general  Ica.-^c-square.s  problem,  it 
is  usually  desirable  to  sot  COLF.PS  at  no  less  than  the  square  root  of 
macliine  precision,  even  if  the  data  are  fully  .accurate.  However,  COLEPS 
may  safely  be  set  to  a smaller  value  if  it  is  known  a priori  that 
th.e  residual  vector  will  be  essentially  zero,  or  if  st-ahi  1 izing  lineai 
constraints  are  to  be  added  to  the  problem.  Even  in  these  cases, 

COLEPS  should  never  be  less  than  machine  precision. 

In  addition  to  the  infonn.ntion  provided  by  COLEPS,  it  seems 
reasonable  to  include  in  the  test  for  "negligible"  some  quantity  b.ised 
on  rounding  errors  that  occur  during  computation  ot  the  reduced  matrix, 
since  clearly  the  calculated  values  will  differ  from  the  exact  values. 
The  backward  error  analysis  given  in  Lawson  and  Hanson  (,197A)  includes 
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the  following  result.  Let  the  k Householder  trims  formations 


Hj  be  applied  to  tlie  m-vector  a,  and  let  a be  the  computed 


transformed  vector;  then 


in  - Hj.  •••  all,  ^ (6m  - 3k  + 40)kE  Hall,  + ^'(g  ) , 


where  C is  the  precision  of  the  floating  point  arithmetic. 

This  bound  is  used  in  HRF.hL  in  the  following  way.  At  the  begin- 
ning of  the  (k  + l)st  stage  of  the  reduction,  the  tolerance 
is  defined  for  the  j-th  remaining  column  as* 


(j)  = (6m  - 3k  + 40)ke  ’ 


wtiere  is  tlie  original  Euclidean  length  of  tlie  full  column. 

If  tlie  length  of  the  computed  j-th  remaining  column  is  less  than  G^(j), 
the  column  is  considered  to  be  "negligible."  Note  that  if  the  computed 
length  exceeds  Gj^(i),  the  remaining  column  could  not  be  zero  with 
exact  arithmetic.  This  test  is  conservative,  in  the  sense  that  it  is 
based  on  an  error  bound  which  is  inherently  an  overestimate;  thus, 
a column  may  not  be  negligible  even  tiiougli  its  norm  is  less  than  the 
given  tolerance,  however,  .as  indicated  in  Section  4.3,  it  is  believed 
that  in  close  cases  a questionable  column  should  be  viewed  as  negli- 
gible, thus  leading  to  a smaller  value  of  the  estimated  rank. 
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Thu  usual  ilufinltion  oj  "ne};l  tgible"  in  IIRT.DI.  involves  a 
tolerance  tliat  is  the  lar^'or  of  COLEPS  am!  the  tolerance  based  on 
backward  error  analysis.  It  is  pos.sible  to  use  only  the  value  of 
COLLPS,  but  this  aiturnativu  is  not  recommended  except  in  extreme 
i-iises  for  the  expert  user. 

Skipping  a t ransforni  ition.  Tf  the  current  reni.iinin,.  column 
to  he  reduced  is  initially  in  the  desired  form,  nr.  Householder  trans- 
formation should  be  applied.  To  te.st  for  tliis  possibility,  the 
Sfjuared  length  of  the  column  is  computed  as  the  squared  length  of  the 
subdiagonal  portion,  plus  the  squared  diagonal  cl.’iiient.  If  the 
suhdlagon.al  portion  i .s  exactly  zero,  or  if  the  computed  squ.ire  root 
of  the  squared  length  ol  tiie  column  is  equal  to  the  tomputed  square 
root  of  r!u>  squared  diagonal,  the  column  i.s  consideied  to  be  alreadv 
in  a suitable  form,  and  the  transformation  is  skipped.  This  is  indi- 
cated by  setting  llV'KCL(k)  to  zero  if  the  k-th  transformation  f-om  the 
left  Is  skipped. 

Recalculation  of  sqiu’ro<^  len^tji!^  of  ^ejnnininj^  columns  . When  a 
column  interchange  strategy  is  used  in  IIRHDL,  at  every  step  the 
'largest"  remaining  column  is  reduced  next.  Anv  nie.isure  o)  the  size 
ol  the  column  will  iiulude  its  length;  wc  now  consider  the  procedure 
for  computing  the  lengih.s  of  all  remaining  columns. 

The  length  of  tlie  remaininu  column  selected  to  be  reduced  must 
always  be  re-computed,  to  maintain  th’  risjuisite  level  of  accuracy  in 
the  computed  t ran? form.i t ions . However,  if  the  lengths  of  all  remaining 


I 
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columns  arc  re-computcd  at  tlic  k-th  staj’i*  of  rcciiirinr.  an  m hy  n 


matrix  (m  > n),  approxlnuitc  ly  (n  - k)(m  - k)  extra  adili  t lons/mul  t i- 
pllcatlons  are  required,  wlucli  will  increase  by  a factor  of  about 
1.5  the  number  of  operations  involved  in  the  total  reduction.  To 
avoid  tills  additional  work,  it  is  possible  in  theory  to  obtain  the 
squared  lengths  of  the  rem.iining  columns  at  any  stage  by  modifying 
their  lengtiis  at  the  previous  stage.  At  the  (k  + list  step,  the  j-th 
remaining  column  is  simply  a t ransforma i i on  of  J-th  remaining  column 
at  the  k-th  step,  without  its  first  component,  i.e.. 


Remaining  column  Transformed  Remaining  column 

at  stage  k remaining  column  at  stage  k + 1 


Bec.iuse  Householder  transformations  preserve  Kuclidean  length,  the 


squarcii  length  of  the  transformed  column  will  be  invariant;  thus,  the 
length  of  the  rem.iining  column  at  the  (k  •»  list  stage  is  given  hy: 


1^11“  = + •••  + = (/.;)“  + •••  + ixjr  - (x;r 


- (xj)"  = H/.ii"'  - (zj)"  , 


I 
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J'UP'iiriNE  tli-'^PL  ( N"  W, 

1 TNTtaC,  SC\LE,  N't-ASK, 
:N'"Fr,i;s  wr.iw,  n „ir. , 

INTEOi';^  IP'iHf-  (NCCL) 
L:ni'>.AL  IM.F'-.: 

PC'JPI.L  Piidl^-'  FC.  N toir’ps 
Doi’bTi,  ri't-ci.-. rc';  OFvivni' 
1 WJF  K ( L) 


SC  [,  NPir,  OPV,  C0LEP5,  r.^OTOL, 
HV  'Cl,  I FiF*',  W-lF'f,  I tpr  Tf  ) 

I,  NFANK,  LFI'FCF 


N:'L),  S'AI,!-'  (NCOI)  , hveckkjoli  , 


c 

c 


c 

c 

c 

c 

c 

c 

c 


THT  JlIjSOLi  'NL  Hhini  flC?  ti  ' :JSEH  CL  DC  F FJCllCII'^V  THE  LFFT) 

r£''inN£.D  tc  c\f?y  j.  ih!  p ''piictic n-  ok  a cenefai  “eal  nfow  '^y 
'iocL  r.A’-.'Mx  Ol'V  T' 1 ijFF.r,  p;  r HZ'": PAL  of  uffes  'rcTANr.iii  ap  fof"' 
'I'FC'IC.H  .,:>t  L I :A'.  It  'J  K r-  E K t'"  E M'" 'i  O’"  OPTHPG'SAl  ( H'"' ME  F liOIP  £ P) 

■"■•A  VSKCPr.  jN  [ 'F".  "Me  FEPnCTl  *’  ^TLT  ''FFM:vaT'='  (A)  WHEN 
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1 EAKT-Sgu  A F SS  SOLUTION  CF  A GFN’='FAL 
Y.INFD,  H MNF 'P-r  ETEP  YISED)  . IF 
"'NT  T''  P"  =»ANK-D'£r  lOIFNI  , THE  AC  OC  .Y  P A N Y IN  G 
ijf'-r  TO  CAFFY  C'J"^  A FUFTHhR  FEP'ICriON  BY 
w^'F(  nv  V- iCNG  ON  IMF  ’’IGH'!-.  I*!  r":S  CASE,  THE 

3''I'I'”T/f;  can  hf  C''''rui "h®'  ''utpu'" 
I'.UT  IHI  SUrF'iiTINt  F:-’'JRT,  to  OBTAIN  AM 
■ r''P'^:jN  OF  CCYFLFIE  '' F"  H ^ G ~ M A 


WT"‘H  '<FFDI  apf  IWC  K''NDK; 


L "C  I'iTERCHANGF  C'lTJFNS  AS  TH'' 

IK  t'A'''^IFD  OUT  (S-''-'  THE  A •''CO«  I- AN  Y T NG 
C.  KO'ISS:  ;N)  . IMF  :nTFPC!!ANGE  STIATEGY 
••FP''CEP  C'^l-TIYN  AS  .HE  N ' v ' 'Ol'IYH  I)  BE 
•j  ri  T?,',  - •'FASH'*  !'  ’ST7"’,  AS  PISC'ISSEO 


: ; ’ -.IFY  ^HE:  HE"  '' 

. g: I A : FI  • : n i 


NCI  ’■N''"FPCMAN  SFS  AF'"  ' 
II  'F  'FALS'=“ 


L IMG 
'0  -'OITI* 
' i.  A '’G  ;■  s : 


YUS"  !F  YFAS'IFEU  IN 
AN"'  .'H-N  .A  t''''LI)rS  I 


'.'S  ".HE  ''LI.^WIkg  IHPE'^  “F1SUFES 


pPEP  T"* 

' K " F>  E SI’’: 


7'. 


T'^S  OPiniNAl  LENGTH 


I 


C jF  A CCHJ/.N; 

J fA)  "tii  FiATlO  jF  IT:;  CGSFF'VT  IEVG7H  T'^ 

C (A  FELaTIVE  AEA  S'lih^  ; 

:T3  C'lhREAL  LEMGT'f  (AN  ADSGLUTE  •EASMSE)  ; 

(C)  TH.-:  ut  ’■"o  CDFI-EK:  IENGTH  I'.  A POSITIVE  SCALTVG 

FAC'TGiv  FFIViDtL'  FV  T"'  USE''. 


I NFGP.I  A 1 ION  CCNCi.RM.;G  fiE  H ,:  :JSEH  C I.  DE  F T F A N ' F'^P  A T ION  S APPLIEO 
FtiO’*  IHE  LtF.  IS  SI  ''FT  -l-LIVi  T!!'  PIAG'^NAI  r?  C'V,  AN''  IN  THE 
AUXll.'A'^y  V-STCi'  tiVSv.T  fFULL  OETAILS  AFE  GTVFN  IN  THE  FXTEFN'AL 
L’GC'IF.FN ’’M’Iv..N)  . IH:.  'tpi  r " F A FI I P T''  WHIG')  IF!  'IPIGINAI  rATFTX 
HAS  HEFN  lEU'JCEi:  IS  .•TGFr'.O  IN  i'li.'  C S F RE  ?P  ' T I NG  PoP'"I''N  F QIV. 
COI.tlFN  IN.  EI.CH  AnGLS  AFf  "'IC  FPEO  IN  THE  A F I Y TPEPK. 


TIiF  F'PEAi,  FAFA'-ETti.'  ^F  ■.  F VFC"  AFE: 


'NTEIEP,  INPUT  "NLY. 

:he  .<i)'iJ2F  ;f  pi.  ws  f ^!>e  ''atetx  o?v. 


rusi  E''  .GT.  E, 


NCCL  - 


INTEGER,  INPUT  wNLY. 

"Hr  .Ji’M’EF  .,F  C LUYNS  iF  T'!'-  EATEiy  OFV.  ‘'UST  HE  . G T . P. 


‘•DIF  - 


INT‘-:GF.  F,  INPU-'  '^NLY. 

■’Ph  JEl^-AFL'P  i.:  » L'I  -,t-NS:ON  :f  the 
SUl'P.UOGFAr'.  ''ilST  El  .GE.  'IPGW. 


'•A'^'^ix  CFV  "N  TH’='  ''ALIING 


('  F V - 


CC'IbLE  rHLC::>Ii..  ANFAY,  0 r-  CU'.'ClPTUAI  PTYFNSICN  N PC  W ?Y  NCOI., 
ANP  PECIAiFEU  ):KfN'':CN  ND'Y  cY  NC'"'!;  INPUT  AS?  '"lJ''PUr. 

N INPUl',  C'IV  , ? 'PIX  TO  3E  FZniJCIE,  TC.  UPPFfi  TPAPEEA'IDA 

Kcr-M.  -ai:,  pe;  :flu  upper  T»ArETc:r  is  st'^pup  in  the 

TCPI-ESP'-NDI.-.G  I ft:  N F 0^ 'J  . CO'tF.NENT?  K*1  THFOUGM  NFCW 
F TiiL  V£c,.r  A.3S  C:A''rD  W 1 1 ii  THE  K-T!l  H i P S F HO  LP  E ? TPAV"^- 
='CP.'iA:ioN  AP’^  SIjFE:  FELDW  ''•he  DI<:'NAL  in  the  Y-7H  CCLI"'N 
OnV.  NC.  f h;t  ,:il  '■''lu‘':.'s  of  t’c-'  '".<a p'^e-tp  hay  ooftespov 
TP  A p u,  •»  i<  N F :ri  j ,lm''ns  ■*. f -hf  .■)'’!g:nat,  '•a'^'^ix. 


(T  I 


r-OUt:Or,  PttEC  S , 

(.OI.EPi  Is  A VAl.ll'' 
"H”  -I ex:  cx;.u'in  *’:■ 


c 

r 


Cl?  J'JSS  I:,  N 
'SA:  e EH'ULl) 


.irn"  only. 

i'f  :j •••;;)  ir;  oEisFr  sing  fc,HEr'F;p  -he  eite 
■E  I EC  'NEGi:;"HF'  (see  TH'-' 


" r'  V 


ip"  I ' 

• 1 ' c 


A P..F:nt-'I'N  *F  •5I7F').  GREAT 


CiHCS^NG  rvjIFIE,  ANP  TH''  il?  ' i; 


I M 


o r,  ( , r,  ''  f~  r ri  r>  c f-  n < ■ n n , r,  n n n r.  c,  n r.  o n n <■5  n ri  ci  o o r,  ii  -)  o n r n r,  n n o 


r 1 

"FCEa  TO  FI  \D  THE  EXTlFNAI.  COCUfPNTATTPN  CN  THIS  SffBJ’J’PT  VERY 
^AfiFULLY.  COL..FS  SHul'LD  3E  NO  SPALLFP  THAN  THE  tlNC  E PT  A I K"  Y IN 
THF  DA.'f,  ; TL.ATIVE  T : H F AFFPOPPIATF  SCALING  FACTOR  AS 
MPNTIONEL)  :hPcF  'XOOTOL'.  CoLEPS  SM'^'ILC  FE  NH  S^ALL^R  THAN 

'“ACiiiNt  pPEJisK'.  £,rc:-:rr  ’’v  fxtfe'ie  ci?cm''st ances . if  htedl  is 

OSilJ  IN  Ct  NJU  M : I.'N  WITH  AN  HNCCN  S TR  A T N F r L F AST-SQU  A S E3 
rFr;,LF«,  THIN  IT  I"  FrFC'FNTLY  WORTHWHILE  TO  USE  A CONSERVATIVE 
VAT  fi  F.R  C!L£PS  (^.G.,  TH'’  SOHAPE  ftOCT  GF  MACHINE  PRECISION) 

Ti.  AVLiO  NU'IrlFICAI  Ox  FF  I CML  1 1 ES  , 

rcnrcL  - 

I EGER,  INPUT  < NT  Y. 

EOPIjI.  is  an  IN^i-GEP  ELAG 
VAIUE  jF  '■.vjJTOL  I'.TICATFS 
FACH  CC  L'lTN;  (2)  IHE 

'•r’SLlGIi.LF'. 

IF  FjniJL  = 1 F -1,  'HE  SIZF  CF  A PEEAINING  COLUMN  IF 
LFFINLb  AS  .HE  FATl-’  F T^S  CtI=r.'"''T  LFEITH  TO  THE  ORIGINAL 
LF‘'.I’H  Or  THE  F'JLI  COL"r.N. 

IF  .TJDIOL  2 uF  THE  r;p  ft  REGAINING  COLU-'N  IS 

DEFINED  AS  ITS  CUt<FENr  TENGTH. 

IF  T.iD.  >L  = J . n -H,  TI■•L  SIZ''  "'F  A PFIEAINING  COIUHN  IS 
DIFINL'''  AS  THE  PA’lj  ■ V IT'S  CUFFENT  lENGTH  I"  THE  SCALING 
FAC.JF  Fh.VxDEL  ^'{  IMr.  HS'^'=  F,F  T"AT  CCLU''N. 

FCSinVE  VALUEr>  ' F v-jL-ci.  ’^NDICAT''  THAT  THE  TEST  F'^F 
' NFGLIGLxL..  ' IS  L'.'.FL'r  G .V  THE  '' F JU'  DSip  - npo  v IDFO 

''CT.E.HA  NCr.  CjLIPL,  *Nr  A QU'VTITY  H'SED  'N  THE  BACKWARD  EFF^” 

AVALYS:>  ^F  H . 'i:  I li(  I.r..F  P A N 3 EOF  " A T TO  N £ A VECTOR  fSE'' 

lAWoJf.  AND  HAN  I xAST -SCH  at  ES  ff:~ELE''S,  P*’ cN-:- 1"' E - 

HAIL,  l;7H,  CHI'T:'  T,  " EIAILS).  N'GATIVE  VALH’^s  CF 

::.o::at..  that  cmy  :t:.''pF  i?  to  ?f  "'■ft  af  a test  '"'ir 

'»'F  I [ ' . ’■"  X'.D'.gl  leg;  'TV?,  COIEU  EIIST  sp;  VO'’- 

AFGA'IVi:. 


WITH  TWC  PUFFCSES:  (1)  -^HL  ABSOLUTE 
THE  “EAStPE  CF  'GIF’"'  TO  BE  USED  FOP 
''F  «''DTOI  INCILATEE  THE  DE'^INITION  OF 


' *’T  0 y'"  - 

lOGICAl,  IN  PI  I NLY. 

IF  :NT--<  13  li'i;  , .li"’  COLII"!;  If  T'EFCMAN  IE  ST'’ATFGY  DESCFIBEO 

AP  VE  -..Li  '•  ULlO.  'I  T'I''TRC  IS  FALFF,  TFF  COIHI-NS  WILT  B F. 

FFDUOF.  .1.  ’■'Hr  siv  N .trF'-. 
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AVfAY  WT!.! 
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FhO.iCEJ,  .•,;.L  VoO  C S ~HF  'JPIOIKAL  COl'IFN  ''PD'^PTvn. 

UT-F  - 

ro'Jb.L  pi-t  ’ijr'’  v:  ••'■  r,  'f  i?ng:’«  nol'i.  'jfjTP'tT  o»Ji.y. 

'HE  ij  ''u  F")"  : F' ;;i- c:  AT  f .^'''fage  rtiPis'; 

F’<’CJ  i .S  r (.'■■VF-''-  I'P'*''  ::  f'H 'HMI  >X-,  G’.V’S  FT 

leng:.>i  •.■  .h  h-"'!  :'lu''N  ''-at  wva  PEDijc-En. 


LEPrCr  - 

TN''E;tt,  G'lTP'i 
A N’  '••!  y i<  F I.A , 


MY, 


5^•^HL£.  VAL’IFF.  AFZ  «5  ’^OLl'WS. 
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I^BPJF  = C : NOFKAL  T FF ” I NAT  ION , NC  FFPCPS. 

lEPBOh  = 1 : tlTHFP  N = CVi  O'’  NCOl.  TS  N V -PO  S T T I VE , OR  r.ODTOL 
IS  MiiGrtllVE  AND  COLEPS  IS  NEGATIVE. 

lEFPOF  = ? : 'lOD'IvL  ^ .1  OP  AND  ONE  OF  TPE  ELFKENTS  OF  THE 

SCALE  ARRAY  10  N C N- P’  5ITI  V E . 

lEBROR  = J : INTER:  IS  AND  ATL  COLUMNS  OF  THE  TNPtlT 

rATRIX  ARE  IDENl  IS.'.LLY  ZFRO. 

lEPROR  = U : INIEPC  To  FALSc,,  AND  AT  SC'IE  STAGE  THF  NEXT  C'AHIHN 
TC  oE  RECncEb  ■ N FGl.  I G I"  T F • . 

lEFFOr  > 10  : A*'  rSPGR  CCCMP'-FD  DOPING  EXFCOTION  OF  PMHLTC. 

TI!IS  SliJULD  NiV'’F  IIAPPFN.  I'"  IT  DC3S,  CHECK  THE  FPFOF  IN  H1ULTC 
C )PREh  P .NDI  NG  . TH:  VALUE  (LFRFOP-10). 


••  AtJTHJio:  .AFGAPHI  li.  UPIGH'’,  SIFVEN  C.  GLASSMAK 

oYSIwES  . PTi:*  IZAT  ION  LAPCFATCFY 
!-EPAR''r.LNT  ry  TSFAIION'S  '’ESFAPCF 
IAN'’.  FD  UMVEPSITY 
S i A N F^  R : , CALIF'-FMA  94  ICS 


■**  DATE: 


rSCEMJFF  197  7 


rtCIAFATIIN  -I-  LjCII  VAFIADL-S. 

IKTEGL.-*  I,  IfcFF,  ;S«VF,  J,  K,  KPI,  LEN,  -AXCOL,  FN-.IN, 
DCIiSLi:,  ii-ECIS.'.ON  <.NuFE,  Cr'riCL,  EPSrCH,  tFFCT,  P"* 

PIEIj,  F.iSl,  SAV-'X,  S..AVE,  T'^IPCK,  VLEKSC,  VSTIA'I 

CFCLA.'ALIjH  v7  ■'TAM/APD  FIJNOTTG'NS. 

integer  lAeS 

rOUPT.E  PI.ECISICN  '^SCF" 

■‘CH  :j  TH-;  f.Ln'.IVI  rCISlON  of  the  FIG'i'TNG  point  afith 
CA'^A  LPSHCH  / E.27I'-1g  / 


PTEPS  - i^SwH  (KPSTCiM 
CF.P''  ^ [ .,gR:  ( r^r KPS1 


I 


c 

r . 

c 

c 

c 


c 

c 

c 


I.EFFCS  = 1 


JHECK  FCy  I^:VALID  VALdZS  CF  INP!!T  FAFA^^TFFS. 


IF  (NFjW  .LE.  0 .OE.  SC  U .IF.  0)  FFTtJ'^'! 

IF  fr.OJTCL  .LT.  u .AND.  I'CLFPS  . LT . C.Ol'  + O)  FETf'FN 

FISC  "HC  .'1INI.'!Dv  NCC'H  AND  INiriAlIFF  THF  TPF?*'  ASPAY. 

IF  (NR..  W .NF.  St-.n  '^r  TI  TO 
MNKTN  - N't'Ji* 

GO  10  20 
1o  'iSF.IN  = 

20  CC  NT  I SUE 

DO  30  T = 1,  SI  CL 
IVZM.iL)  - I 
3 0 C-'MINJO 
N r A N «■  = 0 


C 

c- 


:RF  C-FIGINAL  -ICJArED  0 L'Jr.N  LFNI'.IHS  COTF'Jtp''  ISO  STORED  IN  THE 

'IVECI.  MI'AY.  THE  IFIMSAL  lENIIHS  AFT  S"'''FEC  IN  ' H '='  WORK  AipfAY,  »N'' 
lEI  AINED  :Hr-vUGhDU:  J '.i";  1 0 N ov  hfEDL,  to  PS  rSFC  IK  ;:ONJ'IN~T10N 

jn»  ’ESIC  I.N  VoL N ; i-AIKVapc  Or-FIF  ANAIYS".S  »^''DNTr. 


r 

li  o 

J = 

1,  NCjl, 

CA 

U it  I 

N S 0 ( N F ' . F 

, 0 

PV ( 1, J) 

, iiV 

FCL  (0)  ) 

P.N  (J) 

= DFOF- 

(TV 

^CL  (J)  ) 

CO 

CCNTI 

ter: 

'*0 

"L 

= ! Ac 

0 (K  J'l  :CL) 

or 

(5), 

70,  iOT 

«!!■ 

N 

•• 

t.  J 1. 

1 - 1 

0 

H y.  ^ hi 

INAL 

COLUKN  IF’JGTFS  S'^PV”  AS  CAT.ING 

? AC 

- ^ 

n c ^ 

SO 

cc 

MI 

li  J h 

»>  /> 
t ■ 

h '' 

.)  - 

1 , i ML 

.OCALE(JT  = V^rK 

(OT 

n n n o r,  o n o 


WHEN  Mcn.  )L  = 2 F -2,  AN  ABSOL!I'I>=-  CFTT’^FICN  CF  517’='  IS  USED,  AND 
I HE  S>.  Al.INii  FACS -PS  iJET  10  'JN'IIY. 


70  CCSTINJE 

nc  P 0 j = 1 , 

JCOL 

SCALE  (J)  = 

nO  CCNTTNUL 

1 . 0 P + ■ 

r,~  TO  110 

■'  WH^’N  '".Cr'lJL  ==  1 /P  ■!»  USE^'-SllPPLIFC  fCALE 

C IF  ANY  SJ-iLE(.l)  I^  N-t.-P  SITIV^,  AN  FFPCF  '’^11 
C 


cACTors  A&f.  :JSEn. 
IS  TAKEN. 


9 0 CC‘'TIS  IF 

I^rpprK  = I 

LC  KO  J = 1,  NC-L 

IF  (SCAi-L  (■))  .IF.  C.Cn  + 0)  PE"'UFN 
1G9  CC'^INIL 

1 10  CCNTINIE 


NA^'CCL  = 1 


IF  N'^  I N I E pen  \ N CES  'FE  Z'^  PF  CARPI  ED  ON  I , 
‘JEL'nCiJ  FIRST  . 


THE  FIRST  CCLHNN  IS 


( . N JT 


C 


NTIRTI 


131 


C 

c IF  IMFP'-UAN.'.-S  .'.I’t 
C '’AYI»''IF  CAIlO  liiNCT 
C 

C 


v CAPPTED  OUT,  CETEFIINE  "HE  rCLl)’'N  I''' 


RE  AX  = 

0 . " D ♦ T 

OC  1?0 

0=1, 

!.JCl 

FT., 

J . ■=  W' 

,K  (.’)  /S  • !F  (01 

1 V 

(F . 'oT 

. L = . P.'A  1 0.0 
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I 


* 


C LC'^P  JIATK^-cN:  k i is  the  rMN  Lcor  c;ep  cchjens, 

-o  tE'-riCii  .li:  '-Al-rx  T-  'IPr'?  TFIOV^.ULAF  FC  F ‘i  PY  APPITCATI-N 
C H'.UFLHCLOiti  IhAF  SFOrF  A"  ' JV.’^  FSCF  ""HE  LEF'^. 

C 


n r,  o •*  n n n o 


I 


I 


I 


c 

L 

c 


c 

c 

c 

c 

.!  20 
C 

c 

c 

c 

I." 

c 

c 


c 


r 

r 


C 

C 

c 


COr.l'Ui'''  *H’  SC’IAFEr  LlNr.TH  OF  THF  CDLmN  TT  PF  FKDUCFO  TV  TWO 
PArTb/wiTI  A ;^LTA?A-’F  CA  t.CUT.ATIOV  ' THr  5 U'^  D I A GO  N A I T-OfiTTOV. 

VSLilAG  = C.0D*-C 

IF  (K  .EO.  NS  .'il  ';C  T’.  220 

I ?’.  = NHOW  - K 

ca:.u  l£uso(ftn«  gpv(K  + i,K),  vsoi’.g^ 

F'lit.HEH  PF.TcCII'V  AGAINST  UNDERFLOW  TAN  PE  IVCLGDED  H^PE  IF 
K E«.  JSoAH  Y, 

\ 

VI_.iSv  = VGUlAG  ♦ v'PV(K,K)**2 

■X-iZ  .FETHEF  IHF  rOL'I^^'  3R  r-RSlCSF  IS  • N '^GLI G TP  LE  ' , tlSTVG 

T»E\rI  :L?  ''il  _fF.  I?  -'0  "Y  THE  VAEII’’  OF  VorTOI,.  ^ 


THE  EPOjE  J'  J.<1  B’.oEP  ON  ’’HE  RACKWAFF 
I'Hl'  .ilGINAL  -■  L'jH.N  .■:iGT‘’S,  WIilCH  Af 
A ? r A Y . 


SPFOP  analysis  P,  EO'II  FES 
STOFEC  IN  THF  W'^  P K 


I Ff  N ■ ( K-  t)  * (*' * NP  ~ ♦ UO) 

EEFJ':  = lEF. 

T)LJ-n  - :.PFTT*EFi;:o"*WOrK  IK) 

COLiOL  = 0 ;i.E  T‘' • ST  ALE  (K  ) 


Itorci  .GT. 

"t:?  JN  3>U.NO  AS 


, GjE  TM"  IAFGCP  CF  ''TLiOL  SN  C THF  3AJKWAPD 
A lOLEFANCS.  OTH'^'P^ISF,  'I^F  ONLY  COLTOL. 


I ■=•  (ECDTlL  .01.  . A NP. 

CN'jrvI.  = Dr  JP.I  ( VLt‘’SU) 

I*^  (CN.F''  .LE.  :0!’CL) 


•'LBCK  .GT 
, T*^  3 1 P 


C0LT''L)  C^ITOL  = '’’CLOTK 


I 


THE  K-rii  C >LU''.N  IF  Nu"  ' S'-CLIGIGLF'.  INCFErpN”  the  fank,  an*^ 
pYoTFLJ  WI-'I  li--'  [ FOMCTIv^N  . 


SFA.NK  = N-  A NK  ♦ ) 

!'  VEJL  ( <)  =0  .'Jl  » ■) 

I '^2  I'i"  K-0'  ■’P'.‘..forfa:  t:  V tan  lf.  skifped. 

IF  1 ^ . t'..  . N ?1-  /,  ) G c " ' ) V C 

(Y  Jlc'-.  . -0.  '•  T I 2 «c 


82 


L'M_  HJ'UL'l-l  K, 
IF  (LFtrO'  .NE. 


K , NH 


NCCL,  CSV,  CNc?'^,  HVFCL(K), 


T ' :20 


2ij 


C 

C 


C^N’’  N I " 

1 ( K , . v . ;* - L ) 

IF  (iSlIrC?  tiC  I 


3 00 

2‘^  ' 


!'•’  N'l-  r. 3CHAVO^. 
PPti'Ai'.  A lo  \ F.'P  I 'i 


Al  l'.  T'l  Pf  CAP’’ IFF  ■''MT,  S’'T  i*;;.  XICL  TO  K+ 1 IN 
M >T  Tlf*?  THFCnr;;i  THE  L'T'’. 


riX^  CL  ^ K ♦ I 
rr  . 0 200 


C'^NIi  N'I£ 


P'-APF  TH"  SOOA^EF  LENGTHS 
r'F  IwM.a:!;:.  inp.-'  ''f  the 
F.FG.’H  . 


;F  TI!£  W;,**;  jvING  CCL'J''NS,  and  obtain 
ppp'icpp  :':i!!'N  OF  vAyiHor  ■:al'=’o 


FP 1 ' F ♦ I 

H':A  ' - I.OD^O 

DO  .3  5C  J = KPI,  L 

■Iv'EJLf.Jl  - HVr"L(Ji  - OFV  (K  , J) 

IF  !Hv-;i(Ji  .1  •:.  o.''o^oi  go  t:  pf-i 

.•^TEo*  = L’Gj?  T fHVF  ■ L un /3C  ALP  (J) 

IF  C TlSI  .L!.  Hr.’.F)  'I  TO  2'^'' 

; .•: AX  - i : 

.'lAFC^'L  - 0 

•■;  ‘N.'.s 


Tl‘.-.  SCAi.Fi' 

lf.s;t'!--  .l 


r '.  il  £ I 

x A :•■ : r ' * 


.NO'  !is  s:roi,n  , fy  cotpaptn-; 

v.TK  TH^  Aii  AL-.G'>'IS'  VAlHr  -"H"  LAST  t:-*-’  TH” 
E;  -FT'  3GPATCM. 


IF  (F  FAX  . >1  • 

L K N ■ • Ti  >3  w * K 

? r A / = 0 . c ♦ 3 


C ^ L 


• VT'^C'  03  T'' 


nr,  c n o n t . 1 ^ r.  ( . r-  n n 


r 1 


rc  2b0  J = KE1,  NCOL 

JALL  LF'ESCC-LN*  OFV^KFI,,!),  HVEl'L(J)) 

p:t::  = lso?  " f hv^’ci  ktm/scalf(J) 

xF  (HTto'T  .Lj,.  PKAX)  rC  TO  ?60 

a/.  AX  = aiEST 

•lAXCuL  = J 

2 to 

3;\V-,  F.X  = F.'IAX 
C 


C 3TArF,''''Nr  JuO  EiS)S  'ili;  rMN  LODP  OVEP  THE  CClUFKS 
C 


c 

3C0  CriTTNIJL 

c 


„.)  "Er-E,  AT!.  ■'.N.TIN  COLUfTNS  HAVE  BEEN  CCN3TDEFED  TO  BE  'NON-rEEO' 


T V p O p = Q 
pc-  IIP 


;.T  STAf'^TONT  21'.',  AT  r SltG^  OF  THE  P^'-DtlCTICN  FPCH  THE  L-FT,  THE 
riP'^T  cCI.UTN  I : 0 TtlSIOLP-'C  'N  EGLiriBLE*  . 


no  CCHTTN  IE 
I FTP''.-  = 4 

IF  (.Nrr.  :nt;fc)  a.TiUPN 
T ^ F F 0 ? - 0 
PF"tIPN 

IT  ETA  TE :,T  J2d,  Ll.FFOF  , F . B ON  FFTtIPN  FhZ"  H^tJETr. 
C 

12B  LEFP'.'P  ^ Ltr,'  -.i'  ♦ n 
FF'-  n ; 

c I- 
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fi.5.  Subroutine  HKKDR 


i 

4 

n 

i 

f 

> 

f 

f 

I 

I 


9 

! 


f).5.1.  Purpoae 

iiie  subroutine  HREDR  (for  Householder  reduet  ton  troni  the  £ight) 
reiiuces  an  r by  n (n  > r)  upper  trapezoidal  m.itrlx  to  an  r by 
r upper  trianylc  followed  by  a block  of  zeros  by  applying  Householdor 
transforniat  ions  on  the  right.  This  technique  i.s  used  to  compute  the 
minimum- lengt  li  least-squares  solution  and  to  form  the  complete  ortho- 
gonal factorization. 

6.5.2.  Description  of  method 
St'e  Section  3.2. 

6.5.3.  Keywords 

Houseliolder  transformatii^n;  orthogonal  transformation:  complete 
orthogon.il  factorization;  minimum- length  least-squares  solution. 

6.5.4.  Source  language 

Fortran.  The  code  in  HRKDR  l.as  been  cheeked  by  tiie  PFORT 
verifiers,  and  is  WATFIV-compatible.  All  v.iriables  and  functions 
are  explicitly  declared. 

6.5.5.  Specification  and  parameters 
See  accomp.iny  i ng  listing. 


1 

I 


i 


i 


6.5. (5.  Error  inciioators 


See  accompanying  listing  (the  de.scriptlon  of  the  parameter 

l. ERROR)  . 

6.5.7.  Auxiliary  routines 

HREDR  calls  the  standard  functions  DABS  and  DSQRT. 

6.5.8.  Program  size 

51  Fortran  source  statements. 

6.5.9.  Array  storage 

There  are  no  locally  declared  arrays. 

6.5.10.  Timing 

The  number  of  operations  to  carry  out  the  ri'duction  includes 
2 

approximately  (n  - r)r  addit ions/mul t ipl ications,  r square  roots, 
and  r(n  - r)  divisions. 

6.5.11.  Further  comments 

The  Householder  transformations  applied  on  the  right  are 
represented  by  vectors  stored  in  a manner  analogous  to  that  described 
in  Section  6. A. 11;  components  r + 1 through  n of  the  transfor- 

m. ition  that  reduces  tlie  j-th  row  to  apprt'priate  form  are  stored  in 
those  sami'  components  of  row  and  the  j-tli  compi'nent  is  stored 
separ.itely  (in  MVECR(j)). 
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w 


I 

\ 

\ 

\ 


I 

! 


■ 

I 


V 


t 


Ik 


As  in  HREOL,  the  Householder  vectors  are  normalized.  In  tliis  case  so 
that  the  magnitude  of  the  j-tli  component  is  the  scaling  factor  for 
the  transformation.  The  j-tb  tr.in.s  format  ion  is  skipped  if  the  i-th 
row  is  already  in  a suitable  form,  based  on  tests  as  in  flu-  reduction 
from  the  left  (in  this  case,  HVF,(lR(j)  is  set  to  zero). 


1 


c 

c 

r 

c 

o 

1. 

c 

c 

c 

V 

c 

f' 

c 

c 

f' 

c 

c 

c 

c 

c 

0 


c 

c 

f • 

c 

0 

(' 

c 

c 


c 

I • 


c 


SUP!-.'(I.  l:^E  HFc;nP  ( N’FANK,  SJTL,  CFV,  VSIP"A,  LFPROR  ) 

IVT?nER  NFANK,  NCoL,  NEIE,  I’^FPCP 

COUPLE  PHECI^flCN  CPVfNDIP,  Nr,^Ll  , VA  1 1’ HA  ( NC  CL) 


;HE  -SUlPjUT'NF  HF7D.S  (F'‘P  H OU  SFH  OL  LtP  r’^n(jrTT''V  p?0.1  IH"" 

Rir.l)-^)  13  L ioIGN^n  ' EDUCr  AN  H^’PEP  T PA  PE  7.  TD  AL  ''A'rP':X  TO  HPPER 
'■NIAVOULAK  FCS?1  'Y  A?)I'CA";rv  Ot  ape  ROFRI  at- LY  CltrSHN  HOUFSHrLDFP 
TRAN'- FCHr.^r  .r  NS  -t.  rnr  -roHc. 

:H'=‘  INrlli'  .''ArRlY,  R V , IS  ASSUflE''  TO  CCNTATN  as  NPANK  PY  NCOL 
t;pr-rp  TFAPlZwiPAL  '’A'lF'X  IN  IT;'  FIRST  KFANK  PCWS.  IN  A RACKWAPD 
-WEFP  TI'FIJIJGH  Dii  ?rwj,  : = NFAN«f  STEP  -1  HMIL  1,  TH’=‘  T-fH 
'lOIIFFMOL  nSF  :«AN  ■ FOr  .'•.A-'i,«l  *0  N ST  F I!  CT  E O TO  ANNIHILATE  EIF.^’^NTS 

NRANF  + 1 .JM  C I W I,  Af^EF  TM’'  {!,')  DIAGONAL  FI.EN’-N'^, 

AFP  L’^AVri  SLEZc-MS  1 T!'.'  U'-.H  1-1  AND  1*1  THPCtJGH  N'^ANK  UNALTFFFD. 
’^ACh  SUCH  T RA  N-f  ' P YAI'I  .A  IS  'THEN  APFLI'='r  TO  N..i«S  1 "HPj'lGH  1-1. 

THIS  rF'}CSDi!.-t  i:  rLcC'irrn  iv  o'^tATL'’  detail  in  l»'kSon  and  hanson, 

S'LV’VG  LTA3T-St  lA  FEE  FF'LLISS,  fP  ENT  IC  r- HA  I 1 , 1<57U. 

HFKDn  -ILL  NORAIAIIY  HI  USED  IN  CON  OUL  CT  ID  N VilTH  MFEDL,  TO 
OPTAIN  "ric,  Cj-:EL2Tc  AL  EACTOF  : 7 ati  cn  a senfpal  real 

’'A’^F.iy,  and  TCCJFEUTF  TM-  M vi-'uy-LFNr.i'r!  IFA.'T-SC'IAFES  SOLUTION. 


IMF  F'F'-AL  IAl^A'■iTFFS  • ^ HPIDP  A'’’: 

SPANK  

.N.ECEF,  INPUT  .MY. 

NFA.NK  13  lU.  ••'f  LV  US  LY  DSrf  = 'IINOD  RANK  OP  THE  INPUT  7ATPIX 
.)FV.  MAN--  GIVES  'HL  NUNBEF  op  p A AS  F CP  * A 'I' I C N S APPLIED 
.A  IHt  ’’IGir,  '..r>  THE  SI’E  OF  THE  '^'NM  ti  F PE  c TSIAN'.ULAP 
.A  A . R I >. . 

•'.us:  Ft  .GL.  1. 

NCGL 

IN-SGI.F,  INPUT  uNLY. 

NC..L  IS  IH-  C''LU'‘N  PIXrNSION  CF  TP*'  PATFiy  QFV. 

.'lus:  fL  .G..  NPANK. 

NPir. 

IN  . EGF.}  , IM  UT  ; Y. 

.NDI*  13  TMF  E>  : 'FNALLY  DECLAPEP  P'W  IIFENsun  OF  IHF 
OFV.  FU.^r  HE  .GF.  \PAN>f. 

CFV 

oc,u.iLC  irr  iF-.'s  .r/.'Piy,  of  dfciaf'^o  dt-a'^nsicn  voix  by  nc^l. 

IMJ.  A.\|'  l'IT^UI-. 
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1 


I 


i 


I 


c 

c 

V. 

c 

c 

r 

c 


c 

c 


c 

c 


V. 

r 


( 

C 

n 

C 


L 

C 

c- 


u 


r 
( ■ 


ON  IM-n:,  0?  V IS  ASS'TM.FD  TO  CCKTSIN  AV  L'FPE5  T!^ A P E ZO  IPA  L 
'.ATi'IX  :n  its  K^ANK  »SD  AIL  S'.' :■< -DI  AOOK  A.L 

.iL-fldt.i ->•  ArK  I.'I.-'OPr''.  ''N  .UITPUT,  THI  DF.STRED  nPFStJ 

?B  lANC:  LA  •<  \A~r:y  TS  ST  fed  in  'IPPE'^  '"OIAN'JLE  ?e  qrv, 

ANL  JL..  ;FANK*'  TH‘>.'»Ua!i  NC  C I CF  7HE  I-T'l  HOnSEHOLDEP 

VL'Ii.f  FT-. '^Zl  IN  CCin-'N?  *iCANK‘1  N ZC  L '' F FPW  I. 


VALPHA 

DCiJCLS  PFFCISION  AEF.AY,  W DIKF.SSICN  NFANK. 

OUTtUT  M.y  . 

YAiPHAC)  .;lk:a:*s  zh"  i-th  o.y'on ent  cf  tpe  vectof  that 

DtllN-'i  THJ  F!  'l.ZiiCr'^eP  T8  AN' SFOP^F  A TIOS  THAT  FEPHCES  FCF’ 

I I -■  TH^  ^^.-'r^T-ZZ  F'  P'*. 


LEFFi'F 


IN  . FGEF  , OUTP'i:  ONLY  . 

ZiSFOF  IS  IHF  FLAG  WITH  U ’’C  SSI  ELF  VAL'.'FG. 


lEFFlF  = C ; NCr.LAL  PET'JFN,  NO  EEFOFF  '•ZCND. 

LUFF  )F  = 1 ; ILLEGAL  INPUT  PAP.A*'ETFF. 

„EF.F  >F  = ? ; '■.  FT  '^lAGINAI  A NO  ALL  ^L'''*FNIS  ""O  IFF  FIGH"" 
AF’"  2"PC  IN  CNF  F "HE  FIPS’'’  NPANK  POHS. 


»**  .AUTFloFG; 


yj.  T 


h ! J *■-  !• 


!'  . u ' I T !' 


3TFVEN  r.  GLASS  NAN 


; IS  r.-’:-:.'’  pti.  itmion  l ••  .Ff'EA  t ■'p  y 

I'll  AFZ'.N^N'  01'  OPS  FA  nr  V S "FPEAPCH 
. i .■  N'  F.  r ■'  'IN  I V f 3 T T Y 
jiANr'i",  ci::r'>pN:A  uuns 


.««  ta^L: 


L'-X'r'''PFF  IP  17 


r'E'-.-LAFATlCN  LOCAL  VAFIABLEf. 


r K 

FW 

T ? 

X , * 

I r 

: 1 , j , 

K,  N 

r »:ir 

PI 

p j'rn  T 

E V 

FECIF 

:g  n* 

I 

V r 

f ■ *- 

r*  ^ 

0 

ESI'.  N,  ''*’'.'1?., 

' 

TFC!  A 

r AT 

J.C  N i. 

? 

r ' NDA  '’P 

FU'J : 

NS. 

L 

w 1 

PEcIt 

•. 

L-AP.>, 

r 

r 

ty:  i. 

t’3  It 

ILL*. 

L 

V 1 LlJE 

- C 

n n 

T r AY  '’TF'-S. 

1 t. ' F 

t;  r 

1 

I*-  (f.ANh  .IE.  C)  FE’UP'. 
nc  SGO  . ^ 1,  N'^ANI 


li 
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o o < ! ■~i  fi  r <~i  o n n n '■)  o n n o o ''i  r:  o r o o 


V»i  LPHA  (I)  = 

CONTTNUL 
IF  {VPANK  .01. 
lEFPOP  = 0 
IF  (VPANK  .tV. 


C . OD*  r 


NDir  . )a.  NPAVK  .GT.  VCOL) 


ncl:.i  feiupn 


PFTWFH 


:nF  L^'ip  sta*h:k:;n~  ^ caffifs  'mt  a ssrpwApn  sweep  over  the 

/.'FANK  Vr' 43  . I = WH.ANK  FT'-'P  -1  iPITIL  1. 


NPNKP1  = NPASri  ♦ 1 

lEFrrp,  = 2 

n'l  ‘=tO  :i  - 1,  KFAVK 
: = NEANK  - IT  ♦ 1 


C'PMPIjrZ  THE 
EIE-i’-’M  ANJ 
COSPGIED 
PE  oK  IPPS:. . 


LENGTH  '■  F TMr  VECTCP  C'^VSISTING  ('.?  THE  CIAGONAL 
EIEK.NTF  VFAWK»1  THPOHGH  FCCl  CF  ROW  T.  IT  IS 
'■W  , PA  FIS  ' CHECK  WHE"’H’^F  THE  TPA  NSF'ip.i  ATI  ON  CAF 


SOM  = C.  CLO 
O'"  3 10  J = S E N K r 1 , N C(  L 

SUM  = o'lM  ♦ Oi  V (1 , J1  **2 
= 10  C''NT1MJE 

I'  (SUK  ..  0.''P*  ) GO  "0  GoG 

Sll;.  ^ SUE  ♦ CNV 
PNOSM  = DSOFJ(FirM 

T^sr  WMCTIl.F  FMOFM  IS  7'"’0.  IF  SO,  FFT'MPN  WITH  LFFrC?  = 2. 
IF  (iNJF.T.  .^Q.  0.2'  f'l  ’^SrUFN 

TPS  : iS5T  NCN-Z.  R,  FLIF'n:  CF  TH^  HCH'-EHCIOEF  VECT'^P  IS 
FLCMtST  I,  WHICH  IS  ^T<•rF^  IS  VALPHAl^l. 

FI  EM  = OF  V (I  , I) 

~SI.'.  N = 1 . Ij  *C 

IF  (Ei-M  .IT.  C.or*Pi  Fsis:;  = -■'•s'^gn 
VALrilAdl  CL-d/PN:?'  » ^SIGN 

.H:-  ..ANSF.FIID  O’^ASriAL  FIErdKT  HAS  FA;M’'U[)E  FS^FM,  «IIil 
IGN  JPP'...)I.E  T .'(AT  ' t THE  'IN  T*’ A‘»  5 fi'P  T F.F  DTAGGNAL  rL'^'*ENr. 

CFV  C , 1)  = - 'SIC;.»FK  IF  '• 

VVJlul  AL!  .jILTENTS  •.  F TH  - MOHS  E HOI  F ‘F  VFCIC^  PY  RSur'', 

C S.  ThA’I  fH^  APSCl'iri:  VLIE  OF  "Id  I-TH  .■''•'P'SENT  "S  ..I'lAL 

C T'  Hii  Nl  H Ihl.I/lN  G FAT"”’  FOF  TH”  T-T”  T F AN  f”OF*',  ATI'^N . 
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I 

! 


c 


520 


r 

! 


c 

c 

c 

c 

c 

c 

c 


c 

r' 

c 

c 

r- 

c 


5.r' 


r>'  5 2 0 K - Nr  N KPT,  Ni,' C L 
JhV  (I,  K)  = tPVC, 

Co.'.::  ..ir 


''POLY  rilr!  Tf<A.N:FCP"r.T:  ON  "’O  TH£  FCW.'  A3CV’!:  r^'W  r. 


ir  (I  .i-c.  n :c  t:  '.50 

i."  1 ^ - i 

E'^  55  0 K = 1 , ^ ’ 

HF  liCO'^EHOLnE.''  VEJT'F 
:n  is  p. 7:n'^  appliet. 
V’^CIOP  ’■3  fPcnTPD 
I e;»a  (11  . 


Tin  Dv  1 PF-.  rocT  (PoiPFn  cf  7 

3C"  !o  •PlC'i  THF  TP  A'!3 ' MI 
;r  FIE.I:".?  p T!'M  H 50  5 7«CL^|  2 F 


wOneut: 

H !.  D : H 
r h h f r- 
SE  rA.'A". Y PECIU.SE  : 


STCFEE  IN  V! 


OTTPRn  --  VA  IPHA  C ) 

■n  5 30  .1  = > K 1 , 

L'CIEFT  = DCT^FC 
JC  N.  1 


* OFV  (K,I1 
fiC'’  L 

* V (1,0) 


OFV  {h,0) 


i.' 

C 


APriY  r.i  n .'I  '^rANSF^'n  rATi  iN  ic  T!!”  fow.  this  is 

o'J  u.VAl  ■ :(  ;•  " I '5  :F  ACTIN';  A "'.l!  Ill  PI:  (GA‘'*'A)  0*^  Ti."-' 

.i.tlSEH  : Ei:  THE  F''W.  ..NCE  ACA’-N,  T«t  FIRST 

JErFr-N,  '•  ’ ” TCESIHC:  EFK  VECTOF  IS  IHEATSD  S b P A F ATP  I.Y  . 


C 

lA.'-C'T  ^ ;)CIir.5  / Df  fS  ( VALPHA  r I) ) 

2l'V(K,i)  - O-RNfS,”)  - ':AE’1A*VALPH)  (I) 

j>  'rUC  .T  = Nt*.  Kpl  , NC  •! 

ykV(K,J)  = CFV(F,J)  - GAr.YA  * 0‘^V(:,0) 
5 'JO  'CO.NTINU: 

r* 

f'  C'.T  Cv  t r<  WC  1 "IIFO'IGH  1-1. 

C 

5 0 O'" :.  T : N •)  i 

c 

FNE  "'i-  L-'.E  CJ2F  Tlir  I .■  A f S F ' H « > T IC  SS  . 

r 

55  0 CONirs’lE 

TPtF'N  - G 
F"':'IF  .N 

^sr 


I 


F] 


I 


6.6.  SubroutJ^ne_  LEN.SQ 

6.6.1 . Purpose 

The  subroutine  l.F.NSQ  computes  the  squared  Euclidean  length 
of  a vector  in  a manner  designed  to  avoid  underflow. 

6.6.2.  Method 

The  component  of  maximum  magnitude  in  the  vector  is  located 

in  an  initial  scan.  Thereafter,  the  squared  lengtli  of  a "normalized" 

vector  is  computed,  where  each  component  of  the  original  vector  is 

divided  by  the  maximum;  the  squared  length  of  the  normalized  vector 

must  lie  between  unity  and  the  number  of  components. 

An  additional  protection  against  underflow  is  obtained  by  not 

Including  a ('omponent  of  the  normalized  vector  in  tlie  calculation 

of  the  squared  length  if  its  ratio  to  the  maximum  is  less  than  a small 
-20 

tolerance  (10  ) . 

The  proceflure  followed  in  LKNSQ  is  fairly  simple,  and  is  not 
the  most  efficient  possible;  a thorough  treatment  of  this  seemingly 
simple  but  surprisingly  complex  problem  is  given  in  Blue  (1978). 

6.6.3.  Keywords 

Euclidean  li'ugth;  squared  length. 

6.6.4.  .Source  l.inguage 

Fortran.  The  code  in  l.ENSQ  has  been  checked  by  the  PKORT 
verifier,  and  is  WATKl  V-comp.a  t i bl  i' . All  v.iriables  and  functions  are 
exp  1 i c i t 1 V declared . 
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6.6.5.  Specification  and  parameters 
See  accompanying  listing. 

6.6.6.  Kr'or  indicators 

See  accompanying  listing.  If  the  squared  length  of  the  vertoi 
is  t<!0  large  or  small  to  be  represented  on  the  machine,  overflow  or 
underflow  will  occur. 

b.6.7.  Auxiliary  routines 

LKNSO  calls  the  stand.trd  function  DABS. 

6.6.8.  Trogram  size 

22  Fortran  source  statements. 

6.6.9.  Array  storage 

No  locally  declared  arrays. 

6.6.10.  Timing 

The  computation  of  the  squared  Euclidean  length  for  a vector 
of  length  n requires,  in  general,  2n  Ci)mparisons,  n divisions, 
n additions,  and  n + 2 mul t i pi i cnl ions . 

6.6.11.  I’urtlier  comments 
None. 
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NATIONAL  BUREAU  OF  STANDARDS 
mcRocorr  resolution  test  chant 


I 

I 


VLEUSC  ' 


SUFFruriNE  LEN-Sy  ( LFN,  VEC, 

IN™EGEH  LEN 

itECiSlON  Vi;C(LLN),  VLENSy 

c 

n 

^ 

r 

C THE  SUBFu'iriNE  LiiSJC  (.'C.  PHTFS  THE  SQUARED  niTl  IDEAS  LENGTH  CF 
C A VECTOr,  UbING  iL^E  SAFEGlIfFDS  '^C  PEEVENT  HSCFPFL^W.  IN 
C [■AETICUIA.?,  iHt  ELE^RMI  F FAXirilH  NCDUUJS  IN  IHF  VECTOR  IE 
C r.'CAr"D  DHFING  AN  IHIi^AL  ECAN.  IF  THIS  FIRMEST  IF  NCN-ZEPO, 

C :HF  SyjAFED  lENGi’H  vF  A *.0?rALTZFn  VF'''rOR  (CO  NET  F DC  TED  PY  DIVIDING 
C EACH  ELEO^NI  -F  .'Hi  OEIGINAT  VFCTOp  HY  THE  .lAXINDI)  l£  COMPUTED, 
r T'HE  SC'JAlvED  LENGID  OF  THE  N IRMA  LI  ZED  VECTOR  ‘"IST  l^S  BETWEEN  UNI^’Y 
C ANC  "’FF  NDiri.iSR  Cr  JCMr:HEM?.  THF  iyTlAprp  lojjGTH  THE 

C OPTGINAl  VECrOE  IS  THLii  GIV'^N  3Y  THE  POOCUCT  OF  IHF  SOU.ARR  ')F 
C "HF  "AXIMD.''  ELEhFHT  AND  I'H..  .SODARF^^  L’=’NGTH  THI  VOEFALIZtO  VECTOR. 
C 
C 

C THE  FIR^.AL  PARAh-TEFF  LENSQ  APE: 

C 

C LEN 

C INTEGEh,  IMUr  1VI,Y. 

C rt!L  LENGTH  OF  T"E  VEC^OF. 

C 

C VL'. 

C FLuATING  F’^TN:  ARIAY,  OF  LENGTH  'LFN',  -^NPUT  ONLY. 

THE  VECIOF  NillS-J  SyUAPED  LENGTH  IS  TO  PF  COMPUTED. 

C 

C VI.FNOy 

C FLOATIN  ; POINT,  CUTry?  ONLY. 

C IME  CC'PDIt:  lyUAFFD  LENGTH  THE  ORIGINAL  VECTOR. 

C 

C 

C 

r 

c **♦  AUTHjPS:  KAFGAREI  H.  WFIGHT,  STEVEN  C.  GIAS.'^MAN 
C 6{Sit''.S  L'FTIMZATION  LAPORATC-RY 

C DOrAFTEFN”  CF  OrERATIONS  RESFARCF 

C SiMF  r ['  'IMVFRS’’TY 

C SIANFOrD,  CAITpoFNlA  9U:05 

r 

V. 

c DiC-MOEi  n77 


c 


Cr''L';)T'  SC'JAFZO  L?NC;Tii  1.  F TtlF  NO  PK  A L 17  F F VECTOR.  » F”STHSF 

FAFCnifA^'P  AGAIN.S;  n *1  L'F  ? FI  ..'K  IS  TMCLUDED  EY  TFSTING  H !I'='T 
T'MF  PA"!'^  JF  ANY  TC  IFF  .TAVir'’^!  IF  VERY  S*’AIt,,  TV 

WHICH  C?''-*’  ir  I"  NOT  iNCl’JDFf  IN'  THF  5’)*!. 


VLFNSO  = D.0D  + ') 

IF  (V'-AX  ,j.g.  '(.OD  + O)  F»"IJPN 
’’V'Y  = 10  L*  Vll/iY 
DO  :i  : = 1,  1-N 

;F  (D,\  c-  (VF.  • ( I)  ' -IF.  TV«X)  GO  20 
rr  VL^(i>/VFAX 

VL’^N.iO  = VI. -NSC  ♦ RATIC^RMIC 

20  CON'^IN:! 

VTFNSC  = V ■!  AX*  VA  .\>:<  V !.'■  V?0 
p r T rj  r ;; 

FND 
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6.7.  Subroutine  MNLNLS 


6.7.1.  Purpose 

The  subroutine  MNLNLS  Is  designed  to  compute  the  minimum-length 
least-squares  solution  for  a given  matrix  and  a single  riglit-liand 
side.  It  also  returns  the  residual  vector,  tlie  transformed  right-hand 
side,  and  .in  estimate  of  the  rank. 

6.7.2.  Description  of  mettiod 

The  solution  procedure  is  described  in  Section  3.3.  MNLNLS 
simply  calls  the  appropriate  subroutines  to  carry  out  each  of  the 
necessary  computations. 

6.7.3.  Keywords 

I. inear  least-squares;  linear  equations;  minimum- length  least- 
squares  solution;  overdetermined  linear  system;  underdetermined  linear 
system. 

6.7.4.  Source  language 

Fortran.  Tiie  code  in  MNLNLS  lias  been  checked  by  the  PFORT 
verifier,  and  is  WATFTV-compatible . All  variables  and  functions  are 
OLplicitiy  declared. 

6.7.5.  Specification  and  parameters 
Sei-  accompanying  listing. 
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6.7.6.  Krror  indicators 

See  accompanying  listing  (cIil.  description  of  the  parameter 
LERROR) . 

^ 6.7.7.  Auxiliary  routine.^ 

MNLNLS  calls  the  subroutines  HRKDL,  IIRKDR,  and  QRVSLV,  whicii 
I 

in  turn  call  HMULTO,  LRNSQ,  lIMULVC,  VML’LVC,  rRSl.V,  and  UNSCRM. 

6.7.8.  Program  size 

16  Fortran  source  statements. 


6.7.9.  Array  storage 

No  internally  declared  arrays. 

6.7.19.  Timing 

The  number  of  arithmetic  operations  is  the  sum  of  the  operatijns 
required  by  HREDL.  HREDR,  and  QRVSLV. 

6.7.11.  Furtlier  comments 

A column  interchange  strategy  is  used  during  the  triangular 
reduction  from  the  left,  and  the-  "size"  of  a remaining  column  is 
considered  to  be  the  ratio  of  its  curri'iit  length  to  its  original  length. 

If  the  least-squares  problem  is  to  be  sol  veil  for  one  matrix 
and  several  right-hand  sides,  the  user  should: 


I 
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i (1)  call  MNLNI.S  only  once,  for  the  first  right-hand  side.  Following 

execution  of  MNLNLS,  the  matrix  A will  have  been  transformed, 
and  information  about  the  transformations  will  be  stored  In  A, 

; and  in  the  vectors  ITEMP,  TEMP2,  and  TEMP3; 

I (2)  for  each  subsequent  right-hand  side,  say  the  vector  C,  call  the 

subroutine  QRVSLV  as  follows: 

QRVSLV(NROW,  XCOL,  N’MAX,  NDIM,  A,  .TRUE.,  .TRUE.,  NRANK,  ITEMP, 
TKMP2,  TKMP3,  C,  XC , RESC,  TERROR), 

where  XC  will  contain  the  solution  and  RESC  will  contain  the 
residual  vector  corresponding  to  C.  All  tlie  parameters  through 
TEMP3  must  be  exactly  as  given  above,  so  that  the  correct  infor- 
mation is  available  in  QRVSLV. 
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I 


k 


i 


f 


I 


S'’JPPOUri  rNLNLi  ( NPCW,  NC3L,  N fl  A X , NDIB,  P,  CClli'Z.  , 

1 NFA:^K  , X,  ITS*;?,  TEflPI,  rEf^P2,  tempi,  LFPPOP  ) 

c 

IMTEGiiS  NFCW,  NCCL,  NMAX,  N D IM  , NPAN'K,  lERFCF 
:NTEGE«  ITE  MP  ( NCCL) 

DOUBLE  PttZw'lSION  CO  IF  PS 

DOUBLE  I‘Si,CIS10:J  A ( N D IX  ,N  CO  I ) , 3(NPOtf),  X(N’COL),  TEMPI  (NHAX), 
1 IE  M£2  (NCC;,)  , TEKPMNCOL) 

C 

C 

c — 

c 

C 'T'Hc  3U3POUTINE  XKLXLS  CCMP'ITES  THE  M IN  TM  DM- I.E  KGTH  l.EAST- 

C SQUAI-ES  SOLUTION  OE  A IINEAF  SYSTEM,  AM  CAN  THFIEFopsj  ?L 

C USED  CO  SOLVE  NO  N - SiN  GUT  AP  , 0 VC?  - DSTE  P M NED  , AND  UNDPS- 

C DETEPMINEO  LINEAR  SYSTEM. 

C 

THE  PPOdLEM  :0  be  SCLySE  BY  XNLNL'^  IS:  GIVEN  A FFAL 
C NSOW  BY  NJOL  MATSIX  A , AN''  AN  NFCW-VECTCR  3 , FTNT  THE 

C NCCL-VCCTOR  X OF  KINIMUM  FfCLICEAN  LENGTH  SUCH  THAT  THE 

C EUCLiriAiV  NORM  OF  CHE  FESIDDAL  VECTOR  (A  * X - 3)  IS  A MINIMUM. 


C 

C 

C 

C 

C 

c 

c 

c 

c 

r 

C 

C 

r" 

c 

c 

c 

r* 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

r 


IN  rrOEF  lu  SOLVE  THIS  PFOP.  LEM,  MNLNLS  CALLS  A SFI 
OF  SUEPCUriNES  CHAT  FECl.iCF  THE  MATRIX  A TC  UI-FES 
TPIANGU1A3  rOSM,  AND  IH^'V  SOLVE  A TRANSFORMED  PFCBIEM. 

FULl  DETAILS  O’'  THE  SOLUIION  PROCEDURE  APE  GIVEN  IN  THF 
EXTERNAL  OOCUMEMA  "ICN  . THE  FOUTINES  FEQUTP'’D  EY  MNLNLS  AS:" 
HMULTC,  HRIOi,  HHEOR,  lENSO,  QMULVC,  OFVSLV,  TFSIV,  '^KSCFM, 
AND  VMULVC. 


THE  FCFMAL  rASAMSTCFS  C'  MNlNi.S  AFE: 


N'RCV  -- 

INTEGEF  , IN  PUT  ONLY. 

THE  NUMESF  OF  POVS  IN  THE  MATRIX  A.  MUST  EE  . GT . 0. 

NCCL  -- 

IN'CiSER,  INCUT  ONLY. 

THF  NUM3EP  O'"  COLUMNS  IN  THE  MATRIX  A.  MUST  PE  . GT . 0. 
NMA.y  -- 

JNTEJEF  , INFU  " CNI Y . 

Tr""  MAXI/ifJ.';  OF  (NPOR,  NCOL)  . THIS  VAIf)E  IS  NFFDED  FCF 
FORTRAN  TO  MI.OP  CHE  CCFR"'''"  DYNAMIC  DIMENSION  OF  TH^ 

tempjraey  arrny  tempi. 

NDIM  -- 

IN'^EGEF,  INPUT  ONLY. 

THF  OLCLARID  SC'-)  IIMENSiON  OF  THE  MATRIX  A.  MUSI  °F  .GE. 
NEC-. 


( 


d 


C A -- 

C DOUBLE  PRECISION  ARRAY,  OF  CONCEPTUAL  OIKENSICN  NPOW  BY 

C NCCL,  AND  DECLARED  DIMENSION  ND:*:  FY  NCOL.  INPUT  AND  OUTPUT. 

C CN  INPUT,  THE  ARRAY  A SPCULD  C'NIAIN  THE  r.ATRIX  INVOLVED 

C IN  THE  LEALT-StUAPES  PF03LFM.  CN  OUTPUT,  THE  MATRIX  A HAS 

C PEES  IRANSFDFMFD,  AND  CONTAINS  IN"CFKATICN  ABOUT  THE  REDUCED 

C FORK  AND  THE  TRANSFORMATIONS  USED  IN  THE  RF.CUCTION.  FOP 

C DETAILS,  THE  USER  SHCULD  REFER  I C EXTERNAL  PCCUMENTATICN. 

C 

C COLSFS  -- 

C DOUBLE  PRECISION,  INPUT  ONLY. 

C USUALLY,  COLEPS  SHOULD  PE  SET  TO  THE  MAXIMUM  0?  THr.  SQUARE 

C FOOT  OF  MACHINE  PPFCISICN,  AND  THE  FELATIVE  ACCUP.ACY  OF  THE 

C ELEMENTS  IN  THE  '*ATFIX  A AND  THE  VECTOR  F . OTHER  VALUES 

C MAY  BE  USED,  BUT  THE  USFR  '^HOULP  "FAD  COCUMINTATTON 

C FOP  THE  ROUTINE  HR EDI  FFFORE  DOINO  SO. 

C 

C 3 -- 

C DOUBLE  PRECISICN  VECTOR,  OF  LENGTH  N FC P . INPUT  ANU  OUTPUT. 

C :N  input,  3 SHOULD  CONTAIN  THE  FIGHT-HAND  SIDE  OF  THE 

C LEAST-SQUARES  PROBLEM.  CN  OUTPUT,  THE  VECTOR  P WILL 

C CONTAIN  THE  TRANS  FOP  MET  PI';HT-HAND  SIPS,  0*P. 

0 

C NPANK  — 

C INTEGER,  OUTPUT  ONLY. 

C THE  ESTIMATED  NUMERICAL  BAN.f  OF  THE  MATRIX  A. 

C 

C X -- 

C DOUBLE  PPFCISICN  ARRAY,  OF  LENGTH  NCOI.  OUTPUT  Ginly. 

C THE  CCMPOTEO  .MIN  IMU  M- LF  NGTf?  L FAS  I - RQU  A F FS  SCIUTION. 

C 

C ITEMP  -- 

C INTEGER  AFRAY,  CF  LENGTH  NCCL.  On-PUT  ONLY. 

C ITFX.P  -ILL  JONTAIN  A PFfTOPP  CF  T;'"'  COLUMN  INTEPCHANGES 

C CARRIED  GUI  DUPING  THE  '’SIANGULAP  FEDUCTICN.  IT 

C COPKESEONDS  T3  THE  PAPAM.’^TFR  'IEF”M'  IN  THE  FCUTINE 

C HRED„. 

A' 

C TEMPI  — 

C DOUBLE  tPECISICN  ARRAY,  OF  LFNGT”  NMAX.  OUTPUT  ONLY. 

C ON  EXIT,  TD.iPI  WILL  CON^-AIN  THE  RESIDUAL  VECTCP  FOP  THE 

C LEAST-SQUARES  PPOFLHM. 

C 

Q 'ppvcp  — — 

C rOUBLL  PRECISICN  ARRAY,  OF  LENGTH  .NCOL.  OUTPUT  CNL'-'. 

C ON  EXIT,  rFMP2  WILL  CONTAIN  INFCR'.  ATION  ABOUT  THF 

T TPAN-'OrMATI  JN  3 USED  D'JFIN'G  THE  RF.''ncTTCN  FROM  THE  TFFT. 

C TEMP2  COF.RESPONDS  TO  T‘-r  ARRAY  'HVLCL'  IN  THE  ROUTINE 

C HFEDL. 

n 

C T E M 1 - - 

C DOUBLE  FFECrsiON  AFRAY,  OF  LFNGTU  NCOI.  ^UTFUT  CNLY. 

C ^N  EXIT,  TEMPI  -ILL  CONTAIN  INFC'''ATICN  ABOUT  THF 


100 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TPA.VSFCSMATIOf-o  (’?ED  PI  RING  I'HE  ?ErtJCTICK' 

FROM  1H5  TEMPI  COPFESPONDf'  TC  THE  AFPAY  'HVECF' 

IN  THE  SCUTIHM  HFEFE, 

LEFSOR  — 

INTEGER,  OUTOJT  ONLY. 

AN  ERROR  Fl.IG,  ^ IT  H THE  PCLIOWING  FOSSTELE  VALUES. 

LErROF.  ^ 0 : MO  ERRORS,  NOPMAL  TE  FfllN  ATICM . 

lEFROR  > C : ERiO?  IN  t'FEDI  (CHEOK  DOCOMFNATION  OF  HPFDH. 


***  AUTHORS:  MA:;i;AFET  h.  WF.GHT,  STEVEN  C.  GLASSMAM 
SYSTEMS  C FTIMI? ATION  LABOFATOP> 
DEEAFTAENT  OF  CPSFRTIONS  FE5FAF0H 
SIA.VFOPO  'JNIVrrSITY 
STANFCRF.  CALIFORNIA  9430' 

DATE:  MARCH  ^978 


C 

C 

C DECLARATION  OP  LOCAL  VAriABLLS. 

C 

LCGICAL  INTEFT,  MINLFM 
INTEGER  MODTOL 
C 


IFTFPC  .TFUE. 

MCDTCL  = 1 

C 

C HFETL  JA.RRILE  OliF  IFE  iir'JCTTCN  CF  A 1 0 OPPfF  TFAPEZjIDAL 

C ECFM,  AND  ClMPUTEO  THE  ESTIMATED  BANK  CF  A. 

C 

TAIL  HREDl  (NRCW  , NTOL.  NDTM,  A,  CCLEF5,  F.ODTCL,  INTI'C,  TEMPI, 
1 NFANX,  TEM2,  ITEM?,  TEMP’,  LEFFOF) 

IF  fLERFOR  ..ME.  C,  FSTIJPN 
C 

C IF  NECBSSAFY,  MSEDF  FFDUCFS  '"HF  UPFFF  TRAPEZOIC  OF  THE 

C ^PAN.^FOE.MED  MATRIY  A TO  AM  ^’PFEF  TFIA'.’GLE  FPCM  THE  FIGHT. 

C 


IF  fNBANK  .LT.  MOOLI  CALL  «iF ’'CF  (N  F A MK , NCTL,  NTTM,  A, 
1 IFPROF) 

MINLEM  - .TFUI. 


C 

C 

C 

C 

c 

c 


OPVSTV  CCMPIEIES  THE  SCinilON  OF  THE  LK A ST-SO U A ' E S 
RPCELF.’l  £V  .FAN  3FOf. MING  THE  ^IGHT-KAND  SIDE,  STIVTNG  A 
TPIANGL'LAR  F-YST£M,  AND,  IF  SFCESSAFY,  IF  ANSFOFMIMG  THE 
.'’CL'JTION  . 

CALL  OR  kSLV  (N.RJ  V,  NCOL,  NMAX,  KDIM,  A,  INTEPC,  M.IKLEN, 
1 ITEMP,  1E.MP2,  rE/F3,  E,  X,  TEMPI,  lEFRCR) 

FFTIJPN 

END 


TEMF3  , 


NPAKK, 


K)l 


6.8.  S lib  rou 

6.8.1.  Purpose 

The  subroutine  QMl'LVC.  applies  to  a vector  the  special  sequence 
of  r Householder  tranformutions  constructed  by  HKKUL  in  reducing 
an  m by  n matrix  of  rank  r to  upper  trapezoidal  form  from  the 
left.  Its  principal  uses  are  in  transforming  the  right-hand  side  of 
a linear  least-squares  problem,  and  in  back- transforming  the  residual 
vec  tor . 

The  sequence  of  transformations,  which  are  stored  in  compact 
form  (see  Section  6.3.11),  can  be  applied  in  either  forward  or  backward 
order.  bet  the  set  of  transformations  applied  on  tlie  left  be: 

•••  Q , 

where  11,  reduces  the  j-th  column  to  the  appropriate  form. 

can  apply  tiie  transformations  in  forward  order  (multiply  by 

T 

in  reverse  order  (multiply  by  Q ). 

6.8.2.  description  of  method 
I'tie  sequence  of  transformations  is  assumed  to  be  represented 

in  the  compact,  normalized  form  described  in  .Section  6.3.11  and  6. A. 11. 
Tile  transformations  are  applied  in  the  desired  ordi'r,  taking  advantage 
of  the  known  structure  of  the  Houseliolder  vectors  (i.e.,  Hj^  docs  not 
alter  components  I througli  k - I). 


QMl^VC 
Q) , or 
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h . H . 3 . 


Kfvwords 


Householder  rr.ins  forr.iar  Ion  . 

h.H.t.  Source  laiiKun;"' 

Fortran.  The  codo  in  QMllVr  nas  been  t-liecked  hv  the  PF'jRT 
verliiti,  aovi  ^.'ATKIV-conipat  ible.  All  var  i .-.lil  e.s  and  t'.nction;;  are 
exT'liciilv  decl.'red. 


6.8.3.  Speci float  ion  and  narameters 
See  acconpanyin^  listing;. 

6.8.6.  Frror  indioatora 

See  acconpanyinj;  Ust'ue.  (the  liescri  pt  Ion  of  parameter 

I.F.KROR)  . 

6.8.7.  Auxiliary  routines 

(.'MUI.VC  calls  the  .standard  fuiu  tions  lABS  and  DABS. 

6.8.8.  Program  si^e 

3S  1‘orcran  .‘-ource  statements. 

6.8. ').  Array  storage 

No  locally  declared  arr.iys. 


10.3 


6.8.10.  Timing 


The  number  of  arithmetic  operations  required  by  OMUTVC  is  in 

2 

general  of  approximate  order  2mr  - r". 

6.8.11.  Further  comments 

Tlie  vectors  that  define  the  Housetiolder-transformations  must 
be  stored  and  normalized  as  described  in  Sections  6.1.11  and  6.4.11 
If  HVKCL(k)  = 0,  the  k- th  transformation  is  skipped  (the  vector  is 
not  altered  by  it). 
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S'll'F  MJri  ’ C.'l'JLVv'  { N.lOW,  V-'ANK,  C f V , HViL''L,  V'=CIN, 

1 '‘ULIN'C,  vLc^ui.  ) 


iNaE'^.Si!  N?.:w,  nfank,  soi.r,  >[iLiyc,  lfffh 

nCHFLE  P'ih.Cl.UjN  crvcin-,  VE-ANK)  , !>VI-TL  (NF  ANK)  , VFCr\ f SitOW)  , 
1 VFCJUT  (NFlW) 


THE  F"Frc'j:'iNr  (,;'”nvc  aipij^’S  a F’^c^rNCf  of  hpiff fclpf,? 

A NFF-'Fi'i  A"  IONS  r'  a st'^hep  :v  f^'vaf;.  ''f  ? a-; 

CF  n£F.  rtlE  .>.ivniNCE  ■ f 'Ff  '*  ATIc  vs  IS  AiSII''FE  T"*  HaV' 

oFVF.f'A'T'^u  Hi  it.j;  s Hi- or  f.r-’Di,  aft  foi.i  rrrAiir  Apf'ur'rfiF 
PFOCf'UF'Z  gIi/zn  ’’.n  f'^p  Fy-K^NAL  I'oct'i '■*.‘1  A': ::)t  ant,  :y  ri:?’ 

Ci'*' A'F’*  1 S rtT  T.'IF  EGINN.Nr;  jV  HPPDL. 


■'!£  FCP'-Ai.  PAS A;i.t;l'ZFo  .F  (;;pn''vz  ’.?f  -- 


N'‘'CW  - 

IN'^LGEF,  I'(FUr  MY. 

:h'=’  NUviriF?  jf  Ct.  y. r.  .jr>;r r i*j  tmv  ayc  cm'p'jt  vfctofs,  ant 

"Hr  ■ F S’,  w.:  . F CFV. 


.•JRAS'-f  - 

rN"‘EGi;?  , INF'::  cmy. 

':H7.  no  -."tF  F H.  T FANFF  ZP  YATT'N’S  in  PF.  APPEI”"!. 


:n’T[  :e:  . pot  cmy. 

"H'-'  )KC:.AFFE  "OW  ’'TMN’FION  F’  ip.r  ••AT'“IX  C’V.  •■•'IS''  Sr,  .Sv’, 
NPCn. 


DOtis:.,;.  P.PFClZI'iN  AFfAY,  Ar  DECLAPT  r'W  r:F.£vSi:.»;  Vr  t - , WITH 
AT  LEA;:.:  NFAHK  CCL'i'N?.  rVF'JT  ('MY. 

THF  'A  MIX'  'jri  ..'J'lT  C'HiAlv  rilZ  '‘''ITP’IT  OF  TIM  SIM  ri  p PPEPT.. 


HVECL  - 

PCIJPLE  MFAY,  OF  LENG"'fi  K3ANK.  IK'iiT  ^*JLY. 

TIM  YKJ..\>r  'IV^TL  ’MIST  i";  V F c;  r V'"' A 1 ' F sy  Snp.PO'ITI  N E HFPPL. 


VECI’-  - 

P-*:...  . . ; 0*  AFFM,  OF  IFNITH  N IC '«  . TNTU"'  OMY. 
"IM  VL_i..P  .V  Hr  ■[  ■■  AM  K S'- ’M. 


•'nilNC  - 

IFM"'.!!,  It.  iOT  TMY. 

• IM  YAlIF  Cr  Ml-;*-  I'riCATFj  IH'='  Cspi’'  is  fc.siCH  THE 
’ = ANoF  -f. '’A  ; • .M  -(•  hi  AFPLIFD.  ' H F VECTOR  '1(11 
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DO  "i  n f,  ' . n n r.  , , I':  n n - i . rj  n n n n rv  c,  o o o o i")  r r , c n r,  n r r.  n r>  n 


CCFnESF^NU:  .iU  TC  -HE  T-TH  T F ^ NSP' OP  r A T I C N F(I^  HAS  7LFr!^  TH 

cc’'Pikc;n:.s  1 rn.r'jnh  :-i,  and  non-7;’F'".3  eisfwheb''.  lft 

C D^NJIF  11’^  Eb  DUCT  MNPAN'K)  * . . . • P (2)  * P ( 1)  . IF 
'II'T.INC  = 1,  IH.  r an:  F 'PI^ATIINS  APF  AmiFr  IN  FCPWAFD  CPD3R, 


?(NPANF)  * ...  * P(2)  * t (1)  • VFCIN  = VECOUT,  OP, 

U * VECIN  - vr.C'f'UT. 

IF  M'JLINC  ^ -1,  THE  I- A NS  FOF AT"!  C N S SfE  AcpiIED  IN  ?HF  EEVFRSF 
rpPEa,  I.F., 


IM1)  * 0(2)  ♦ ...  » F(NraNK)  * VSCIN  = VECt'1",  TP, 
(J  TEANiF-.SE^  • VFCIN'  = VECODT. 


D'nyLE  [bECl'.ICN  A.fPAY, 
IDE  TiiANorCP  l.-.E  VE'.  T F. 
VFC^N  AM  VII'.UT  .'AY  r- 


''F  LENGTH  NFOW.  CUTFUT  ONLY. 

THE  ACTUAL  PAPAHEIEFS  COF FES P JNDI NG  TO 
T H “ S A r E . 


MMCF  - 

IMTEGEF,  OUTFU"  CN!Y. 

A'-  EFPCF  FLAG,  WITH  TH.L  FOLLOWING  PCSSJLLE  VAIUKS. 
loFFCF  - '■  - NC  ''PF't'',  NCPNAL  T E F ^ t N » T I':  N , 
TKSEiF  = 1 - INVALID  INPUT  PA  I A F . 


AUTD.P.):  . AKGAFET  !1 . WIGHT,  STEVFN  i" . GLASSTAS 

SY.STr''.S  ,.PTI.  17  ATION  LAPCFATOPY 
;.cPiFIMEN"'  OF  OPSFATICNS  FLSEAFCH 
•lANF'Fi-  UNIV'PS^TY 
:aNF  FL,  CAMFOPNIA  94  3 OS 

OA'^r;:  F SCNE’' H :-P  1977 


’'FCLAIAM'  N VATTAnTE'-. 

JNIFGiH  .JH,  GF~1,  N'COP 

l)C"r)!E  rFi:Cl..r';N  P'.PF’,  fact 


r,'.  i" 


I 


1 


c 

c 

c 


:a.  , 

rcf-Fi/.  rr-K-i.i.'j 


C ;!‘FCK  F''!-:  tfii.rf  t \ XNTfl  TA  f A rFTF  i-'S , 


C 


1 F.  F P I ' !■  = 1 

IF  (rA;5:j  (.'-IJ  L7 '.C ) . “.  . 1 . OF.  nro,  .t_  c, 

• 'I.  NPAN'K  . f r.  NF  PFl'i?*.' 


. N'’A  NP  , L'C  . ■■* 


" •rT.v'J'v  -P  HA-rVAPP 

'Ft  LI'  A . 'N  ,r  Tifr  71  A •; ' r -)f,  m ^ ^ 

L 


(FIJLT  N 


0) 


]0  TC  1) 


JH  = 1 

PC  T'  tic 

1 0 

J H = N i<  A N !S 

?0 

CCNTII.  JC 

PC  JP  i = 1, 

* i'  . ^ 

V SC  j u 1 { : I 

= V;',c: 

■ ) 

! •> 

ccmtf  r-' 

/ • 

' r r F F F = ,) 

r 

- — - — 

— 

-- 

— 

’ HF 

’OOP  rc.  siAi 

h - ’ 

1 IT 

C 

7C 

P F A P P I r C . 

. '!  ■■  : V 

!_P' 

• PF 

' ' 

A VP 

L 7 ~T)  k:  7 ■'  r,  i. 

C F..NP 

r 

L r . 

c 

(FrrwAPn  APfL.Ct 

: I'  !.) 

X 

FF(  M 

A PF  [ : 'ATIv  N > . 

V7?  7!'E  N'U'^r-FF  Cr  I?AS’?F0P-A;  I-?;3 

’i  ’■'IE  IN^FX  TME  TFA  N EPF-p  V , - J-^,. 

’ '!  ‘^HNP  F »:■■■»'  1 X _ N c A v K ’’p  ’''ll.'  N''  - 1 

N'PA';K  TC  1 IF  »ni.TK’C  = -i  (f.\CKi?APP 


DCTP'rD  = UVECl.  (.TH)  • Vf\  ^tl-"  f JH) 

IF  (JH  .tC.  NPFW)  r.r-  ‘T'o  <^o 

JHPI  - Jf  ♦ 1 

DO  4u  : = .iHF  1,  n:  W 

D.-iiDRO  = DOTIFO  ♦ V E C''IJT  ( I ) •OK  V ( I , 0 H) 

4 0 CON’.'i-.iOt 

■j"  CCNriUUa 

AFPLY  IHf  TRAt.SK'lF.MAiI  CN  PY  SOP"?  AJT I NG  AN  AFirrppiA'^F 

MflLriPl-E  OF  •'hr  V HE  S PC  LH  ? F Vi:cTOF,  WP'^PE  IFF 
COF.PON^M  IC  rEFAUn  f FPA  KA’^FLY  . 

FA.:r  = l)Otfed/uai-..  (Kvi''L  (OH) ) 
v?lju:(jh)  = ii"  (.11  ) - fac-’*'ive:i  (.1H) 

IF  (JM  . £0.  NI<  ■ O or-  TO  7 0 

D'  bO  I = -iSiri,  N'  • 

V.^)iir(I)  ^ v;-  'IIO  - FA-IT-OFV  (I  ,.IH) 

^0  ccn::m)F 

If'  JH  --  jii  ♦ rjT.iv'" 

100  CONI  IN  (If 

I HF  LCrr  'lAj  iefiina  in  ;.oPr»Li.Y. 

I’^EFO.''  = '1 
FFIIM  S 

t NP 


i>.9.  Subroutine  (iRVSLV 


6.9.1.  Purpose 

The  subroutine  QRVSl.V  is  used  to  solve  the  linear  least-squares 
problem  for  a sinp.le  right-liand  side  after  the  matrix  involved  has 
been  reduced  to  an  appropriate  form  through  application  oi  Househdder 
transformations  i)y  subroutines  HREIII.  and  HREDR. 

b.9.2.  Descruption  of  method 
See  Section  3.3. 

6.9.3.  Keywords 

Linear  least-squares;  linear  equations;  overdetermined  linear 
system;  underdetermined  linear  system. 

6.9.4.  Source  language 

Fortran.  The  code  in  QRVSLV  has  been  checked  by  the  PFORT 
verifier,  and  is  WATFTV-compatib le . All  variables  and  functions  are 
ex[)  licit  ly  declared. 

6.9.5.  Specification  and  parameters 
See  accompanying  listing. 

6.9.6.  Error  indicators 

See  accompanying  listing  (tlie  description  of  the  ])arameter 
I.KRROR)  . 
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6.9.7.  Auxiliary  routines 


QRVSLV  calls  tlie  subroutines  (JMULVC,  TRSbV,  UNSCRM,  and  VMULVC. 

6.9.8.  Program  size 

38  Fortran  source  statements 

6.9.9.  Array  storage 

No  locally  declared  arrays. 

6.9.  10 . Timing 

The  number  of  aritlimetlc  operations  is  the  sum  of  the  operations 
required  for  QMl’LVC,  TRSLV,  UNSCRM,  and  VMULVC. 

6.9.11.  Further  comments 
None . 


7-e  , 


S'll-l’Cvl " INii  Ol'V'iLV  ' Ni’CW,  SCOT.,  QPV,  INTPF', 

1 ''’M.-N,  N.-NK,  . HVtCP,  >-■ , x,  if',  lzpoch  ) 

IV'r'mi  MU'W,  LiCl.  N-SAV,  NUT.':,  iJPANK,  I'-^FC? 

I S T t.i  i. F L4.\  . (S i.  '.  T 1 
LCfU^Ai.  Isrs-T,  r.  ^NL  M 

r‘''P'c;  uF.-.ci  f.  '.'PV'Nii*',  ^Jri.)  , pvfci  fNcci) , uvecs  (kcol)  , 

1 AtNjjii,  .'1AX) 


THf  mJjp  'iri.  E i'f  VJf  IS  !'S"U  ?'■;  SOLV.-'  TME  I'PEF?  LfftS'l-- 
SOI'AP"?  "''.rsIEf:,  .>!:n  J-NOf-  (CSV*X  - IT  IC  .■!» 'SU.'lEn 

’’M,'.  r TH?  ^..Vf;->:  JKV  !'A.'  PE-'  'EPLCiT  T^'  T'’ I A N'T  ii  L A ■>  Pf-?M 

f'Y  .n£FLTJ\'  IjU  ^ ■ HOHStH  EL'  I.  F H A'.'.  TFC  ?!•  S TT 3 FRCy,  T»K  LEFT, 

?C£rilI.v  -iT.K  C L'F'.S  IN':  .T-'MEANnE::,,  A fiP  (OPTTTN«LLY)  Tf'A'^  I H' 

"ATS'^X  W.\o  FUill'.iF  I'L'JCir  lY  A r’ F L ICA  T TC'J  C 'r  H C OS  ' H OLD  FE 

:rA  NT  f'."?.';  i : : - m:  " n tiE  ’tgl:  t?  roF  FA-'-Fry  .apptAPFD  '"o  be  mnk- 

PEFTCTE:'!'.  -Hi.  CA.LC''L'.T1''NS  ABE  rAB^IET)  COT  EY  THE  S'!B»o 'IT  INS  S 
B-F'n  A NU  liS-:''’-,  AH.  CJH  F A-'aTLIAPY  hCU'TNFS  (PiSCPICFD  IN  CFTML 

IN  T l‘  ■■  T'U  JO  N . M /,  I ■■  N - F .H^  I ' AN  P '!F  '""’I  . 


c 

I.N  PPUSn  I- 

L V P : H E 

1-; ST-SOUAFF? 

I '?  ■ I E '■  ' ' 

F A f.I  VEN  1 I GHT 

{ ; 

[!Av:  ::r-  vp. -pea 

i*  , i t'  - 

FCI lowing  steps 

*.;iST  E?  c 

MPIFF'  ■'III: 

r 

(A)  AHFLY 

u Pf’J  T 

" AHL  FOF.NATl  ONS 

THAT  'WFFF 

APPLI'Pi  'N 

c 

T'*?  LTFi.  . 

:[ '1  G'l' I 

.0  """  FA'LII’’ 

FPrc'Ci:"  N, 

yiELE"VG 

A TPA\3FC.-r'i  J V-'  'TJf  ;♦!  , 

('  ) "H  LI.-^JAP  SYFTF“  F *y  = .''F,  »HEFF  F IS  THE  N^'ANK 

r-Y  NciAl.r.  UFi  ;-.  1 I A N G'' L .A  F " A T ' * X IN  TB'  UFFF?  LE’^T  COHNTF 
'•  "i.'J  AFA  N;  FCFy.'T  QhV , A N' D EP  COWAIN'S  IlE  FIFST  NFANK 
C''«PJ.Vf  Nl-i  i TH'  VFJTt-P  O'"’.  NFA’JK  I:'  THE  N I""  FF IC' !.  'ANK 

crv,  -.ni.AiEr  o'ifing  the  "apiife  fecoctipn. 


(J)  IF  Ai.  1:  . ’A  N;,F.:FyA"'l'' NT  XEFF  AFFIIED  TC  CFV  ON  THE  FIGHT, 
ArtiY  iiifcs.-  ...  AN.  F''F''A.':t::;:'.  i''  y,  yz?:.rj:ir.  v^y  = x. 

Nj  Tr'A.N  . ).T''A'TrN  VFF-:  AF.r-LIFU  CN  THE  FIOf-T,  X ^ v. 

(P)  IF  t.:-T'!.'..J  FGHA' GI  s XFp'  CA’FII'.P  OMT  COPIWG  "HE 

rET'l.jriv^  t ArV,  1 1. . i r ' F' A "G  ' THF  A F? ' ' F ? I .A  ""  F r .A*' F NT  I 

f r-  V 


I ) jL>x  f j;  i.apn  :hf  of-idua-  y^'otcs  fci  the  "figinai 

rPjULJ.S,  Aii-..Y  c PANGP..  5'-  T-'  Til''  FFSTTHAI  OF  TH"  TPAN"- 
.'•-'F-'E;)  ■.  , APPLY  PtiJ  HPi'IFPH'l  TrF  NS  f OF ’'AT  I GNS 

■fN  FPYF.F;:.  U cl  1 . :Y  . 'N  rpirPICN,  THT  FEET.  HAL  0 TH'' 

"f  A NjECKF  .0  PF^ELF-  PP  7 » P ''  !N  ICyPCN'^N'^S  1 T NFAN’',  AN" 


THF  SLin.MNIN;  CCrPCNiNTc  ATE  EQUAL  1 'i  THE  CCFPESPONDI  NG 
Cr.lPONENiS  OF  0*L*. 


THE  F'^F.IAL  PAF  A ,1  RTF  F r OF  QFVSLV  APE: 

SF'^W  - 

INI’ EG' F,  IN  PIT  :>SLY. 

"HE  N IfIbEa  JF  f QFV,  AND  THE  LFNCTF  OF  THE  VECTOR  9. 

NC-  I - 

INTEGER,  IKPr.  ''NLY. 

TilC  N IMPF’S  OF  C.LUFNE  OF  QFV,  AN  C THE  LENGTH  OF  THE  VECTORS 
IPEFrj,  HVECL,  HVECF  , AND  X. 

NHAX  - 

INTEfFiv,  INPUT  CNLY. 

NiAX  ^ F.  AX  (NF.«,  NC  L).  THE  LENGTH  O’'  TH'  VECTOF  RES. 


INTEGER,  INPUT  CNLY. 

CHS  DCCLAPEC  P''  «i  DIF.SVSICN  ■jf  THE  '•ATFIX  QFV.  HUSI  BE 
.G^.  NFCW. 


DOUBLE  PIECIEICN  APFAY,  ^ DECLARED  'CW  DIMENSION  Nni.H,  AND 
C'NCLFIUAL  SIZE  t.  F . « FY  N''OL.  TNFU'^  CNIY. 

QPV  CGnTMNS  iHt  'RESULTS  CF  CDF.PUTATIOS  OF  HREDL  M'D  HREOR, 
ANO  CvRnESPONDS  To  "H’"  CFTGINAL  EA^RTY  IN  '’'HE  L F AS  T -S  Q'J  A RES 
IFOai  EH  . 

INTFPC  - 

I/GICAL,  INPUT  CMY. 

IF  IN  i-C  IS  .r'U''.,  CULI'IN  INTSFCHA  ENGES  HAY  HAVE  seen 
CARNICD  cur  OUF.SG  THE  FFDUC'’'ION  'PC  A THE  lEFT.  THE  VAI.UF 
OF  IN.il.:  IN  THE  CALI  TO  QFVSLV  S F'' U 1 D B •'  THE  SAHE  AS  THE 
VALUE  ■'F  IN-’-RpG  IN  "HE  EAELIEF  CALL  '’'O  HPETl. 

HIN1''N  - 

T.'^GICAI,  INPUT  CNLY. 

“INi.rjN  SHOULD  bC  GET  TO  .TRUE.  IF  "HE  1 IN  I"  U*"- LENG  T H LEAST- 
''QUAFoS  iCoUIIlN  IS  ’'EE’’'F'’,  AN"'  1 .FALSE.  OTHFFWISE. 

HINL::,N  IS  .IMr.,  :N"E’5C  ''lirULP  HAV'  PSEN  .TFU'. 

N s A ^ K - 

iN:.,Gt?,  I Ni  ui  " v:  V. 

T!!S  LJiI‘'A:‘EJ  RANK  CF  QFV,  AS  CO-FtlTED  PY  THF  FC'iriNE  HFFDL. 
IFF.  FT  - 

’NTSGE?  rtFNAY,  ( F L SGTH  HC'^L. 


o o ''i  r,  n ri  n n . n r'  (<  r>  f ) Ti  o o r>  <■',  n '/  r;  r r.  n n r.  n o 


I 

C TPJi.'l  JoSTAlMi  INPJRr^TrON  ABOUT  IH  i’  COLMKN  INTE9CHANGHS 

..  C .TAY  HAVf  'KT’  fAfFT!'’'  O'JT  BY 

C 

L !ivrc  I - 

C i:''UaL?-.  PPhlJISI'N  "SPAY,  or  length  ncol.  input  only. 

C THE  VETIOh  HVF.L  IS  GENEPATEO  PY  THE  SUBPrUTINE  MRECL,  AND 

C CONOAi-Nj  I NP  JF.r  AT'f  ON  AE'^'l-"  THE  P A NS  F CP  f AT  I C N CN  THC  LEF'^. 

HVEi'P  - 

DOUdLii  PP£,iSI-N  ARFAY,  OF  LENGTH  .S'JCI.  N F 0 T ONLY. 

THE  VCOrOE  HVtGF  IE  GENEt’A'^ED  L Y TH''.  S'lBFCtJTINE  H^ECP,  AND 
CINTAI'!^  I .nF  CPf  ArrON  APCHT  THE  TP  AN  S ' 08  “A  T I C NS  CN  THE 
HI  GUT. 

r-  - 

DOUdLE  t'PEOiSI  N VECT'‘F,  E L FN  G F E N F '' W . INFUT  •.'’0  CJTI’UT. 
ON  EMPY  TJ  OFVSIV,  H SHOULD  CCNTAIV  THE  FTCHT-HAND  SICE 
.'Hi  i.  EAST-GO  I'AP  E.'S  FC^PLEH.  ON  rv-j  wpr.r.  ORVSLV, 

P HAS  A'EN  OVEShHlTTEN  PY  C*P.  WHFOE  C IS  t'HF  FFODUCT  OF 
THE  Ho’k.EH  ...DL  F AN  SFCFH  ATIONS  .’PPIIED  CE  THE  L^FT  . 

Y - 

DOIJOLE  FFEOISTON  ARFAY,  OP  lENGT!!  NC''I.  CUTFUT  ONLY. 

X WILL  i-lN^AIN  (.0''tU''E:D  SOLUTION  CF  THE  L EA  S T-S  O'l  A R ES 

FP03LS''. 

I FS  - 

C0U3LE  PREOISION  ARRAY,  O'’  L''NGTH  OUTPUT  '’'NLY. 

RES  Will.  CCMAIS  THE  rESIPUAL  VECTOR  CF  TH'’  OPIGTN.M. 
I’-'AST-SvUA  FE  S FPC  ^LE-. 

LOFPCR  - 

IN  5GEP,  CUTPUl  ONLY. 

AN  iER.c-  FLAG,  WITH  .HE  FOLLC'WING  POSSIBIF,  VALUES. 
lZi-..iOh  - C - "C  ''TFOPS,  N.'PFAL  T'"’ ("  I N A T I C N . 

LLRPCF  < 10  - '’.RFOP  IW  C^ULVC  (SEE  OrULVC  DOCUE'’ NTAT I ON ) . 

LEFPOF  < 20  - FRF.F  TV  TPSLV  (EE''  TSSLV  COCU  H FN  T.A  T IC  N)  . 

LLrpCF  < ,10  - EPfCF  ’’N  V^ULVC  (iZE  V.MJIVC  DCC  Ur  E NT  ATI  ON  1 . 


C THC  Sr.'BR./UriN''  OPV.SLV  CALI.S  THE  AHXILTAFY  SU  P F Cr  I N’=' S Q''.ULVC, 
C irSLV,  VMIIVC,  AND  iJNSC.'^'*. 

G 

c 

•••  AUTHORi';  UAoGAFrT  H.  WP’SHT,  .STEV'-N  C.  GIASSFAN 
C .SY.-.rEr.S  FTIf' IF  AT’"'!.  LAJ'ipAT'PY 

rf.c>A.»T:iE  S'  'P’^fation'  rfeiarcf 
c SI  \N  K ? r u'.  :v  f .s  * ^ y 

C .,.T\NrrFC,  CAi:  E''F'’TA  Ui4  3''‘ 

c 


1 1 5 


1 


r,  r>  (-J  n r>  f , n r-i  o o o o o n 


rtHrUFAIION  CF  Li-LAL  VAPIAPI^*^. 
INTFf.FH  1,  NFSKPI 


APPLY  w :o  :ii>:  Pir.ir.-'UM'  sirr  vtcrop,  r. 

CALI.  g-lIJ  LVC  W,  N:A?.K,  NPI",  OPV,  HV^^'L,  E,  1,  P,  LFHFOP) 

It  (T.2t(l'OP  .M.  01  S.-irUFV 

FC[vt  :‘iAvc';’Ar  SYsrsfi,  wit!'  th"  ftpei  vp»nk 

corroN^Nis  jf  ri! :;  ii. an  ■>FCP.’'En  e .aj  the  FinH’'-HAKD  fids. 

c 

CALL  r.^i;LV(NEAN  (<,  N'CT,,  » D T *1 , OPV,  E,  -1,  RES,  LFPROP) 

c 

IF  (Uar  C?  . FT.  1^)  TO  Tf  20 

LEFFOh  = L^:pp^.•^  + 10 

CVTIJP;^ 

c 

20  C'NT’NUf 

IF  (.NJl.  INT^PC  ."p.  .NOT.  UNL^N  .OS.  N'"ANK  .EQ.  NCCLI 
1 c",  r . JO 

r- 

APPIY  r'l.'  FAirlX  V T'  IHF  VFCIOP  r FS  , WHICH  CCPTAINS  THE 
C SC  I HI  UN  OF  '.HP  LINtA'-  EYSTE'*. 

CALI,  VIHLVCIN  :C  L,  NP'S'i',  NP^Y,  OPV,  PVFCS,  FOE,  1,  P^S,  LERPOSI 
'F  (IEHFlR  .r,.  '■)  GC  ■^C  JO 

LFFF  OP  ^ I EPI  IP  ♦ 20 
Pri H FN 
JO  CONTI  Nil! 

IF  (INClRC)  OC  T'..  *=0 
C 

C IF  WErh  NJ  I M r PCH  A N 0 , STlP*'  THF  SCnriON  N .Y . 

r- 

rc  'i"  ..  = 1,  .jii 

'<  C)  ^ r Cr  (I) 

'.0  CfK'lNH;. 

CO  t.o 

^ ‘F  in  . if  ::h;,ngk”  wffl:  "jOt,  PF->FrFr.  if’-  ‘?oiiitton  v'-c^op. 

f" 

so  CENT’N'Ut, 

c^il  ii.I'cf.i(i,  n'  I,  insr,  ffs,  y) 
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* 5 


I 


p 

A?  CCNTT.\UE 

r 

oI.FE  .’ri£ 

b 

IP'lAi, 

V"i 

i':0'  “ IHF  TP  ft  ^ 'F'' .7 '' f C F'’OPIF'‘  IN  THE  VFTTO 

C 

F?F.  :Pl' 

I'"'!  t- 

: N ~ 

AN  K 

COfFO'ONTS  THp  I ? ft  m E F^F  ^ F T '^E.Sir'JAL  APE 

C 

AS.prr''..:,^ 

i 1 

E (At. 

:e  y 

7..!^,  ’NP  THE  [>f:  (NFCU-\-h;,nkf  CO“OOH.-NTS 

C 

■\rv 

r V 

ri!'.  : 

•CCT 

fvr Cv -'O A‘,K)  rcrp.. N FST F the  v.sro5r.''o 

f • 

FTCHT-ilANL 

bl 

f',  p 

c 

rc  70  : = 

NPANF 

= 0.0:0 

70  CCNUNJr 


IF  fMiANK  .EC.  NC.W)  r-FCOF*-^ 
NFNK'’I  = tit.ANK  ♦ 1 


^ s V 

K f>  1 , 

NT 

A 

P'T.S 

in  - 

0 

CC' 1 T' 

JE 

AFFTY 

.'hi:  H_i; 

S L b.  w 1 

t Lr 

■ T ■ ' ft  N S 

F'-p'-ftTI  3N 

S I ’■ 

•■»FVEFiE  CPDEP  , 

00  CPTAIM 

THf  P.T 

JO  A! 

^ r j 7 \ 

c 

F H .- 

■’PIGINAI 

F'"  ?! 

f M ^ 

CftT.L  0 

.'iULV^  fV 

3 % , 

S'  ? ft 

i-K,  NC 

OPV, 

Fv"r  I 

, ^00  -1,  FES, 

LEPP-P) 

r-  E ^ r)  P N 

7 VII 


6.10.  Subroutine  TRSI.V 


f 


6.10.  I . Purpo.se 

Subroutine  TRSLV  returns  tlie  solution  of  a non-singular  r 
by  r triangular  linear  system,  where  the  triangular  matrix  is 
assumed  to  be  stored  in  the  first  r columns  of  an  r by  n matrix. 

In  the  current  package,  TRSLV  is  called  by  QRVSLV  in  solving  the 
linear  least-squares  problem. 

6.10.2.  Description  of  method 

If  the  matrix  is  lower  triangular,  fo’o.’ard  elimination  is  used; 
if  the  matrix  is  upper  triangular,  the  solution  is  obiained  bv  back- 
subst itution. 

b . 10 . 3 . Keywords 

Triangular  linear  system;  'I'rwaf'  e .tTiii’  a.  6 : j tituii  ii . 

6 . 1 0 . -'t . Sourct  ' aigu.  f 
Ki'rt  far  . • ■ 

viTi  f ii'i  , lud  i s \ 
dei'  I a l ed  . 

t' . 1 ' I . S . , I . 

s.  I a' 


L 

r 


■ I 

, ! 


km 


I 

I 


0.10.6.  Error  indicators 


\ 

f 


L 


See  accompany i ng  listing  (the  description  of  the  parameter 
LERROR) . 


6.10.7.  Auxiliary  routines 
None . 

O.IO.E.  Program  size 

35  Fortran  source  statements. 

6.10.9.  Array  storage 

No  locally  declared  arrays. 

6.10.10.  Timing 

The  number  of  arithmetic  operations  requiri.d  to  solve  an  t 

2 

by  r triangular  system  is  of  order  r /2. 

6.10.11.  Further  comments 

Because  TRSLV  may  be  cal  led  to  solve  an  upper  trapezoidal 
system  involving  an  r by  n (n  > r)  matrix,  llte  number  of  components 
in  Mie  solution  is  .il  lowed  to  exceed  r;  components  r + 1 through 
n of  tile  solution  are  simply  set  to  zero. 


I 


I 
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Sf)bR0UriN2;  TFSLV  ( 
C 

TN'^Ei^ba  N,  N^JL, 
ncuRib;  PFiEci..:..N 

r 

C 


N,  NCi.  L,  NDIf,  OPV,  F,  INCSl, 

NPin,  IMCSl,  LFRPCO 

OF  V (NT’-*!  , N)  , F (N  ) , Y (NCCH 


LFPpnp  ) 


c 

c TFi;  snp,r,jTiNO  i^Jiv  's  '.isl'l  solve  a tpianouiap  system  of  ltheap 

C EO'-'ATIO'.-'  (‘liHFR  LOU  ; B '?  lipFSP),  WH^'PE  TRIANGfILAP  FATPIV  IS 

C IN  r*'F  '.IPI’EF  LEE’’  N FY  N COPN'^'F  CF  OPV.  TFSLV  TS 
C SrECIFIC«;.I  Y SEiir.NFr  I BE  t'SFP  TN  SCLVIN’G  THE  IINSAF  LEAST- 
(•  SC'IAFFS  LBUBLE.'. 


C 

r 

>• 

r 

C 

C 

c 

r 


C 

r 

c 

c 

c 

c 

r* 

c 

c 

c 

r* 

c 

c 

w 

r 

L 

C 


t: 


:i!E  E',  F’'Au  LAF  AH.-TEF.,  F '^F^'IV  APE  -- 


IN'^E'JEF,  INP'JT  'MY. 

THt  SIZE  ..F  THE  Tl-.:  Al.Gli  I AR  EATRIY. 

VOCL  - 

INTEGEF. , INPtIT  FNIY. 

TH’"  .VUF.HEP  -F  UNKNO-iN',  ANF  -"HE  N'lFHEF  CF  CCLH''NS  IN  OPV. 

FH'':  BE  .CL.  N. 

SDl  “ - 

INT.-GEF,  INpn:  ONLY. 

'■■BE  OECLARFO  P DIEl-N;  ION  OF  Of  V.  PE  .GF.  N. 

CFV  - 

r;n!Li.t  B.'c,c;s:o;  ar'^av,  of  ''Eclapec  ^:'!Ensi''n  not'*  by  ncol. 

:HF,  iFiAoG'iLAF  MF.y  u,tO  TN  TFSLV  I'’  JCNl.MFEO  TM  III^  nPPFP 
TFFT  OCFNEF  CF  CFV. 


DCaLt  FPECISIlN  Vi  'TOF,  F LENGTH  N.  INFUT  CN’LY. 
■^HE  RiOli.-H.\Nn  SITE  F '"HE  '"YSTET  CF  ""  C '’AT  I"*  NS . 


I N~.Ev..EF  , INF'J  ’ (MY. 

"B"  VAL'I."  TN'-.l  I'-DirAT’^S  WHETHFF  B E YAM'IY  I*'  LO'WZF  ''•P 
PP"':'R  'i  . AN  :il  LAP  . IF  T ' G " = 1,  TB""  ''ATPIV  "5  CCMSIOEFED  PC  BE 

!CW:-;.(  . AN  ifi:  AP  , AGP  T I"  ? J A D F L I il  N .A  T I I N 15  ‘"A'^FT^'n  'IJ-". 

IF  xNC.'''l  = -1,  '"Bi  .('..MY  '5  tJPF"""  fFTANGHLAF,  .A  N RAOKWAPP 
.'OI'i^iON  IS  USED. 


PCHBi,t  EFECIMCI  Al-PAY,  OF  lENGTH  CPTfHT  CN'IY. 

'IHF  Flf,..:  t CC'irCN'":‘^„  F V APF  FHF  SOI'ITTCN  '■'F  'YHC 
'•KTANGIJLA?  .->YSiLr  , -.Bn  .HE  '•EHAININS  roirrN'NTS  A^F  SFl'  T'O 
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no  r -•)  r>  o r.  O iT  n O n n 


I 


lEF  T-CT 
IN 
AN 


TEGES,  OUTPUT  CM.Y. 

SRBCP  FLAG,  WITH  THE  F-^LLCWING  r.'!A  SINGS. 

LSRSGF  = 0 - NC  •'fPrES,  NORMAL  T F P “ I N A"  T C N . 

T.'ERRCF  ^ 1 - input  PAFAMSTFF. 

LErFOK  ~ \ D.AGCI-’AT  EIFMfNT  ’'F  THE  TPIANGULAF  MATPIX  IS 

21  • ■. 


AUTHORS 


MAFGARE:  H.  WriGHT,  sieves  c.  classman 

SYSTEMS  •'PT  I--'I2.ArrCN  LAPC-AICRY 
DEL' ASir  LN  ■ OF  SPEP.MICNS  AF.'^FAFCH 
S f A NF-  Pn  MMVFP  5TT  Y 
SIANF:  FL',  CAIIF-'PMA  94  ^CS 


♦**  DA’^E; 


EtCEMPEF  1G77 


PI.''!  AFATI  . N :F  local  VAFIABLES. 

ICTSGES  I,  :i,  I'^lAnT,  K,  NIMI,  NLCOl 
LiC'IELE  LFECISI'^N  yV.aL 

CECLAF.ATISN  of  MANP^''!)  FUNCTICN'S. 

1N~EG2R  1AE.5 


TEST  FCF  JFFOS  IN  INFni  PA r A « ET2 F S . 

I E t F 0 S = 1 

1'=’  (N  .LE.  C .CP.  NCOL  .L".  0 .OR.  N .G", 
IF  (TAjS  (iMCSL)  .N'f  . 1)  »E?UFN 

INITIALIZE  IP.t  Y VICI’.:  •’!'  7IrI. 

nc  io  i = 1,  <coL 

Y (!)  = 0.0'J«-C 

10  C^’TlNJt 


NCOn  FFT'JFS' 


C SFT  ur  IKOIlES  J S.IVF,  ^ITHTR  \ LCWEP  IRTANGIF  (F^HWAPD)  R ' N 
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o n n n n ■;  n ii  n 


1 


UPFFB  TFIANGLE  (bACKWU-D). 


IF  (IN'JSL  .LT.  0)  GO  T.,  20 
IS1APT  ^ 1 
GC  T:'  jo 
20  TElAF-r  = N 
10  CCNTINUi 
K = To  tap: 

LSFROr  = 2 


THE  LCCt-  rC  SlA.i.'lENr  100  TONS  ''VFP  TH'"  'JNKNOhNS.  K GIVES  iH.E 
INDEX  Cr  IHt  UWKNOiiN  CME.>LN:IY  HFING  SCIVET  PC  B . 


nc  ICO  NLCOP  = 1,  N 

IF  (0BV(K,K1  .'^Q.  0.0P*0)  petupv 

yVAL  = t (K) 

IF  (NLC'I  .EC.  11  GO  TC,  ftO 
NCP.  1 = NLoJP  - 1 
I = IS. AFT 
no  50  ::  = i,  mfi 

YVAL  = YVAL  - G?v  fK,  I)  • Y (I) 

: = I ♦ INCFL 
50  C(  NTINUE 
hO  Y (K)  = YVA  L/I.BV  (K  , K) 

K = K ♦ IMS' 

1'-J  J-NIINU., 

0 ■'HF  LOOP  HAS  TELHINAIED  NOF*'ALLY. 

C 

l'=’FFCR  = 0 

FE'^'iJEN 

END 
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I 


fi 


1 

6.11.  Aul>.rou t i ne  _ Wrs JA;  ) 

1 

6.11.1 . Purpose  ' 

I 

Subroutinf  TRTSLV  returns  th<  solution  of  a non-siniular  r x r | 

linear  system  wlierein  the  matrix  is  the  transpose  of  the  trianjiular  ^ 

matrix  stored  in  the  first  r coJi:mns  of  an  r bv  n(n  2 ^ matrix. 

TRTSbV  may  be  useful  in  solving  linear  sy.stems  related  to  tin  orthr.(;oaal  | 

f ac  t or  i za  t ion  of  a j’iven  matrix  — for  e.xample,  in  findinj'.  a t olut  'on 
to  a set  of  r linear  equality  constraints  involving  r variables. 

6.11.2.  Description  of  method  i 

The  standard  methods  of  solving  trianp.ular  system.s  are  used, 

as  in  TRSI.V.  The  only  notable  feature  of  TRTSLV  is  that  the  matrix 
involved  in  the  linear  system  is  the  transpose  of  the  indic.ited 
triangular  form.  Thus,  if  the  original  matrix  is  upper  triangular, 
forward  elimination  is  used  in  solving  its  transpose. 

6.11.1.  Keywords 

Triangular  linear  system:  transposed  linear  system. 


6.11.5.  Specification  and  parameters 
See  accompanying  listing. 

6.11.6.  Error  indicators 

See  accompanying  listing  (the  description  of  the  parameter 
LERKOR) . 

6.11.7.  Auxiliary  routines 
None . 

6.11.8.  Program  size 

37  Fortran  source  statements. 

6.11.9  Array  storage 

No  locally  declared  arrays. 

6.11.10.  Timing 

The  number  of  arltlimetic  operations  required  to  solve  an  r by 

2 

r triangular  system  is  of  order  r /2. 

6.11.11.  Further  comments 
None . 
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o r o T r-  r:  n o '■)  r,  r-;  o i~,  o o r)  o o o o r o o r j ri  o r.  c j r ; ri  T;  o '■)  n 'i  o 


Hf'--  'M  . ; Nt 


: o 


*•  » 


Y,  LfTFCK  ) 


I 


.L  V f 


1 


INTE''.F„<  N,  nd:;*, 
PP  £.  i;-  lu  S 


:ncpj.,  !."''P"ip 

Or  V (VI  ■,)  , L fV  ), 


V (M 


fflF  S'tUPjU'^  INE  TrISI  V IS  L'OfP  f'  SOLVE  R SYSirf*  CF  LTVPRf  FOUATrOSS, 
WH??F  THE  KATFxX  TO  FF  lUrP  TS  TtXNSPC'‘F  'F  A TPIP’.'GULAP  “ATPIX 

(ICKFP  OF  JPFlP),  Sr^ELD  I'.  !»'••  'TPEP  TFKT  K hy  K CCP*’EP  Or  OFV. 
Ti-IFLV  ’S  Sp  tCII- ICA  LxY  OEOIGNEH  '•"'R  TFIFLE*'?  ASSfCIATFO  WIT.! 

AN  CP  H''G.N'AL  F A C T ; I I Z A T Ii.  N ^ OWV,  KHEF''  II.NSA»:  FySlF-^S  .AFISE 
INVCLVI'lj  b.IU  THE  T f 1 AV DL  A F '!«TP1X  ANF  ttc;  TFAPSr''Sr. 

T”"  PSAi.  ihD  CF  "'FT3LV  APF  -- 


TMFOEh,  INPD:  only. 

"■HE  3I7-;  OF  THE  rPTAPGlJLAF  HATRIX. 


INT'EGEP,  IN  FIJI  CNIY. 

"HF  OEOLAPE;-)  POW  D'^lNflON  Of  V . r.UST  pF  . GE.  N, 


ypv  - 

rO'IfwF  PRECISION  AP.'!AY,  OF  CECI.APFr  I' T H:=’N'SIC  N NDIM  PY  N. 

'HF  "P  lANGl’LnP  'firlX  dSEO  IN  TFTSIV  IS  CCNTAINFC  IN  THE  UPPER 
l’='FT  CcFNEP  C F vFV. 


DOnELE  PhECISI'.  V V.J"'?,  ■)  "^  LENGTH  V.  TNIUI  ONLY, 

TflF  kIGHT-!!ANli  ^IOF  o F t['E  SYSTFr  CF'  F C'i  A'T' I C N G . 

INCFL  - 

INTEGER,  INPUT  N L Y. 

"^HE  VhlU;:,  >.?  INCSL  INDICATF.S  WHFTHEP  "fiS  YATPIX  IS  L ' WFF  CP 
"FPER  iRIASGULAF.  IF  INC.SI  = 1,  TF'C  HATH  IX  IS  OCNS:DER!='r  TO 
PF.  LSWEn  Til  ..A  NGU  LAP.  IF  TN^SL  = -1,  THF  HATPIX  IS  AN  UPPF.R 
tPIANGL;:.  'JOIF  ‘"HA?  TPF  '*ATPTX  USED  IN  SCIVINT  T !!  E LINEAR 
^'Y'JxE''  IS  It:.  TfAN  Pfsr  Cr  "HH  TYPE  OF  TriANGLF  INPICRTEP  HY 
"IF  Vrt  L’lE  CF  TNCSI  . 


l''^UHL,r.  PPEi  ICIC;;  PPRAY,  F LENGTH  N.  C’J  T PUT  CNLY.  | 

-RF  FI.tST  N CCHfCN.-.  N'''S  u Y APF.  THF  SOl'iTICN  C’='  TH?  3 

TP’ AN  GULAP  TY.ST'X,  ^ ND  rErAIIHV',  CO'*  P-^N N K-'=‘  .SET  TO  f]l 

c 

C l .'S  PC?  - I 


1 

I 
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TN^EviEFi,  CUIPUT  .TFI,  Y. 

AN  lNF-F  flag,  W::h  "ItE  F)LLOWINr,  FEANINGS. 

LEFFJa  = 0 - NO  EFF'‘PF,  NOF.rAL  TEP  I NA  TIC  N. 

Lc-t  rCn  = 1 - II.  VALin  IJJPUT  PAPAPETFF. 

L’FFC?  = ^ ClAGr.NAL  ELEMENT  OF  THE  '^FTANGULAP  MATPIX  IS 

TfP'''. 


AMTHuRS: 


***  DATE; 


'APGAPET  H.  WPTGHT,  STEVEN  C.  GIASSKAN 
SYSTEMS  F PTIFI7.\TI0N  LAPCPATCRY 
CLPAPTMENT  OF  OPEPATICNS  RESEAPCH 
SIANFOPD  UNIVERSITY 
ir\NF.,Fn,  CALIFORNIA  943CS 

rECEMBF?  1977 


CECLAPATICN  C?  LOCAL  VAFIAFLFS. 

IMEGER  1,  II,  IM,  ISTAPT,  K,  SLMI,  NLOOP 
DGIIPTE  PR  El:  SION  YVAL 

DFCLAPATIJN  .F  SXA1.LAED  ItlN''TIONS. 

INIFGER  lAfcS 


T^ST  FOR  elf:?  in  INEUI  pa:  AP*"'FFS. 
LFFFOa  = 1 

IF  (N  .LL.  0)  PET'J'^N 
IF  (lABSCKCSL)  .NT.  1)  RFTUPN 

INIIIALIZE  IHt  Y VECT'~F  TC  7EP0. 

O'  10  1 = 1 , N 

Y f I ) = C . C L ♦ C 
10  C'.NTINUr 


otl  Ur  INDICES  Tj  solve  ZI.H''R  a LOWFF  I’-IANGL'  (F^PWACO)  JF  AN 
'JILEF  ’•RxANGLE  (PACN./.fD). 
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6.12.  Subroutine  I'NSCRM 


6.12.1.  Purpose 

The  subroutine  UNSCRM  constructs  a re-orclerecl  n-vector  by 
applying  a specified  permutation  ior  its  inverse). 

6.12.2.  Description  of  method 

Let:  tperm(i)  denote  the  1-th  component  of  the  permutation; 

x(i)  denote  tlie  i-th  component  of  the  original  n-vector;  and  y(i) 
denote  the  i-tli  component  of  the  re-ordered  vector.  If  the  permuta- 
tion is  applied,  tlie  vector  y is  defined  by: 

y(iperm(i))  •>  x(i),  i = 1,  2,  ....  n 

If  the  Inverse  permutation  is  applied,  then 

y(i)  x(iporm(i)),  i = 1,  2,  ...,  n . 

6.12.1.  Keywo  rds 

Permutation;  re-ordering. 

6.12.4.  Sourct-  language 

Fortran.  The  code  in  UNSCRM  has  been  checked  by  the  PKORT 
verifier,  and  is  WATKlV-compat il)le.  All  variables  are  explicitly 
declari'd  . 
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b.l2.'S.  Spoc  i f icat.  ion  and  p.irainoti  rs 
Soo  u-comp.iny i !ig  listing 

6.i2.(i.  i rror  indicators 

None.  Iho  iisoi  is  rosponsihlc  lor  tlv  consistoaoy  rf  the  sno!  i 
tied  permutation,  i.o.,  oacli  integer  between  1 end  n sliould  occur 
only  once.  Tlie  actual  parameters  corresponding  to  the  origina’  and 
re-ordered  vectors  must  not  be  the  same. 

b.ld.7.  Auxiliary  routines 
None . 

6.12.8.  Program  size 

18  Fortran  source  statements. 

6.12.9.  Arr.iy  storage 

No  lor.illy  declared  arrays. 

6.12.10.  Timing 

Tbe  re-order  in;’,  requires  n indireit  addri'sses. 


6.12.11.  Further  comments 
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6.13.  Suhroutlno  VMl'LVC 
6.13.1.  Furposo 

'I'lie  .subrontino  VML'LVC  applies  to  a vector  the  sequence  of 
r Householder  transformations  constructed  by  HRKDK  to  reduce  (from 
tlie  right)  an  r by  n upper  trapezoidal  matrix  to  upper  triangular 
® form.  VMULVC  is  used  principally  in  computing  the  minimum-length 

i least-sc|uares  solution. 

j The  sequence  of  transformations  may  be  applied  in  either  forward 

or  backward  order.  Let  V be  the  orthogonal  matrix  which  is  the 
product  of  the  transform.ations  as  constructed: 

I V = H^  . 

I 

I 

t 

where  the  vector  corresponding  to  11^  is  zero  except  in  positions 

i,  and  (r  + 1)  through  n.  \TilH.VC  can  apply  the  transformations  in 

T 

forw.ird  order  (multiply  by  V),  or  in  reverse  order  (multiply  by  V ). 


6.13.2.  Description  of  method 

i'he  transformations  are  assumi'd  to  be  stored  in  compact,  nor- 
malized form,  as  described  in  .Section  6.5.11;  they  are  applied  in 
the  specified  order,  taking  advantage  of  the  special  structure  of  tlie 
transformations.  If  llVF.CR(i)  is  zero,  the  j-th  transformation  is 
sk i pped . 
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6.13.3.  Kcyword.s 

Housetio Icier  reduction  from  tlie  right;  mi nimum-lengtli  least- 
square.'^  .solution. 

6. 13. A.  Sc'urce  language 

Fc'r trail.  'I'lie  code  in  VMld.VC  has  been  checked  by  the  PFORT 
verifier,  and  is  WATFI V-compatih le . All  variables  and  functions  are 
explicitly  declared. 

6.13.3.  Specification  and  parameters 
Sec  accompanying  listing. 

6.13.6.  Frror  indicators 

Se'e  accompanying  listing  (the  description  of  the  parameter 
I.ERROR)  . 

6.13.7.  Aaixiliary  routines 

VMl'lA'C  calls  tile  standard  functions  lAbS  and  1)A3S. 

6.13.8.  Program  size 

16  i'ortran  source  st.itements. 

6.13.0.  Array  s ! o ra ge 

No  internally  dei-Iared  arrays. 
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i 6.13.10.  Timing 

The  number  of  arithmetic  operations  required  is  of  approximate 
I 

» order  2r (n  - r) . 


6.13.11.  Further  comments 

f 

None . 
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